/* File: AquaCompl.c

**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht, the Netherlands and the BioMagResBank, 
**  University of Wisconsin-Madison, USA.

 
   AUTHOR(S) : Jurgen F. Doreleijers based on programs by Ton Rullmann
**  VERSION   : 3.2
**  DATE      : January 2001
   PURPOSE   : Calculates completeness of NOE restraints and can present
		lists of contacts having different properties.
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"
#include "AquaTypes_Assign.h"        /* ADDED; also includes the nrutil stuff */

/****************  function prototypes  ****************/

#include "AquaRange.h"
#include "AquaFuncts.h"
#include "AquaFuncts_dr.h"
#include "AquaStrucset.h"
#include "AquaFiles.h"
#include "MenuFuncts.h"
#include "Functs.h"
#include "AquaHow.h"

#include "AquaCompl.h"
#include "AquaCommon.h"
#include "nrutil.h"

#define DEBUG_INFO true                        /* Printing extra info */
#define NOT_A_NUMBER -99.9

/****************  globals  ****************/

extern Boolean  iaf, batch;                           /* import from Functs */
extern int      terseness;                            /* import from Functs */
extern Range    *range;                               /* import from AquaRange */
char            *program_name = "AquaComplete";       /* export to Functs */
char            *program_vers = aqua_version;         /* export to Functs */


enum {  opt_ext, 
	opt_sho_selonly,                                    /* For Ranges */
	opt_savpseu,  opt_delredu,         opt_delredl,     /* For Redundancy */
	opt_usedist,  opt_thresh, opt_detail, 
	
	opt_atominp, opt_use_intra,  
	opt_numb_shells_prs, opt_min_dist_prs, opt_max_dist_prs,
	opt_numb_shells_rst, opt_min_dist_rst, opt_max_dist_rst, 
	opt_double_count, opt_anal_ssa, opt_avg };

/* BEGIN_EXTH */

int main( int argc, char *argv[] )
{
    char         *fct_name = "AquaCompl main";
    Strucset	theSet, theSet_ref;
    Mol		*molptr;
    Model	*models = NULL;
    DistinfTyp	ftype = inf_extended;
    Qtyp	qtype, qtv;
    int		s_count, count;
    Filnam	dinf_filnam, drst_filnam, dval_filnam;
    Filnam	filnam_s, filnam_i, filnam_r;
    Card	ident;
    int		item, exit_flag;
    DistdataList inf_list, val_list;
    Menu        *oldmenu;
    Boolean     done = false;
    FILE	*fp = NULL;
    char         temp_value[4];
    
/* My own variables */
    int		struc_count;
    int		number_not_observable, number_not_in_range, number_ssa_double;     
    int		number_intra, number_low_only, number_nobs_atom;
    int		dist_count, obs_atom_count;
    int		obs_restr_count, red_restr_count, no_restr_count;
    int         obs_pair_count,  red_pair_count,  no_pair_count;
    int		multi_restr_count, *rev_mapping_table;
    Distinfl	*dist_infl = NULL;
    Distrst	*dist_rst;
    Distvall	*pair_val;
    Distval	overall_completeness_factor, thresh, avg;
    Boolean     *obs_atom, save_pseudos, del_red_l, del_red_u, use_distance;
    DRestraint  *obs_restr, *pair;
    DRList      DR_list;
    char	*atominp, *PathName;
    Filnam	filnam, rng_filnam="-", FileName;
    Extra_Distinfl  *extra_distinfl;
/*    Distinf	*dist_inf;
    Distvall	*dist_vall;
    Distval	*dist_val; */

    enum { item_range, item_do };
    static MenuItem menulist[] =
    {
        { item_EXEC,  "exec",  none, on },
        { item_SET,   "set",   none, on },
        { item_range, "range", none, on },
        { item_do,    "do",    none, on }, 
        { item_QUIT,  "quit",  none, on } 
    };
    static OptionItem optionlist[] =
    {
        { opt_multi,       "multiple",    logical, on }, 
        { opt_batch,       "batch",       logical, on }, 
        { opt_terseness,   "terseness",   integer, on }, 
        { opt_prtaq,       "print",       logical, on }, 
        { opt_prtres,      "printres",    logical, on },
	
        { opt_rest_sel,    "rest-select", oneword, on }, /* For ranges */
        { opt_sho_selonly, "sho-selonly", logical, on }, 
		
        { opt_sterlib,     "stereolib",   string,  on }, 
        { opt_ext,         "useext",      logical, on }, 
        { opt_guess,       "guess",       logical, on },
        { opt_delredl,     "delredl",     logical, on }, 
        { opt_delredu,     "delredu",     logical, on }, 
        { opt_savpseu,     "savpseu",     logical, on }, 
        { opt_usedist,     "usedist",     logical, on }, 
	         
        { opt_use_intra,   "use_intra",   logical, on }, /* General to AquaCompl */ 
	{ opt_detail,      "detail",      logical, on }, 
        { opt_thresh,      "thresh",      real,    on },
        { opt_rstinp,      "restrinp",    oneword, on }, 
	{ opt_prtlib,      "printlib",    logical, on },
        { opt_atominp,     "atominp",     string,  on }, 

        { opt_min_dist_prs,    "min_dist_prs",    real,    on }, 
        { opt_max_dist_prs,    "max_dist_prs",    real,    on }, 
        { opt_numb_shells_prs, "numb_shells_prs", integer, on },
	
        { opt_min_dist_rst,    "min_dist_rst",    real,    on }, 
        { opt_max_dist_rst,    "max_dist_rst",    real,    on },
        { opt_numb_shells_rst, "numb_shells_rst", integer, on }, 

        { opt_avg,             "avg_power",       real,    on },  
        { opt_double_count,    "double_count", logical,    on },  
        { opt_anal_ssa,        "anal_ssa",     logical,    on } 
    };

    static Menu menu =
    {
        "Main menu", "command", "Give command:", 
        sizeof(menulist)/sizeof(MenuItem)
    };
    static Menu options =
    {
        "Options menu", "option", "Give option:", 
        sizeof(optionlist)/sizeof(MenuItem)
    }; 


/* initialize */
    range = NULL;
    inf_list.count = 0;
    val_list.count = 0;

    ProgramInit( argc, argv );
    ProjectInit();
    StrucsetInit( &theSet, "theSet", true );
    StrucsetInit( &theSet_ref, "theSet_ref", true );
    SetMenu( &menu, menulist );
    SetOptions( &options, optionlist );
/*
    SetSubMenu( &options, opt_rest_sel, &menu_rest_select,menulist_rest_select);
*/
    ModOption( opt_rest_sel,	"both" );
    ModOption( opt_sho_selonly, "true" );
    ModOption( opt_multi,	"on" );
    ModOption( opt_rstinp,	"aqua" );
    sprintf( temp_value, "%d", terseness );
    ModOption( opt_terseness, temp_value );    
    
    
    PathName = GetVar( "AQUADATADIR", true );
    sprintf(FileName, "%s/StereoLIB-aqua", PathName);
   
    ModOption( opt_sterlib, FileName );
    
    ModOption( opt_thresh,      "9999.9" );	/* For fixed distance removal */
    ModOption( opt_savpseu,     "false" );
    ModOption( opt_detail,      "true" );
    ModOption( opt_use_intra,   "true"  );

    ModOption( opt_double_count,"false" );
    ModOption( opt_anal_ssa,    "true" );

    ModOption( opt_atominp,     "ob_standard");	/* Observable atoms */

    ModOption( opt_min_dist_prs,    "2.0"   );  /* Distances pairs */
    ModOption( opt_max_dist_prs,    "4.0"   ); 
    ModOption( opt_numb_shells_prs, "4"  );	/* Number of shells used */

    ModOption( opt_min_dist_rst,    "2.0"   );	/* Distances restraints */
    ModOption( opt_max_dist_rst,    "4.0"   );
    ModOption( opt_numb_shells_rst, "4" );	/* Number of shells used */

    ModOption( opt_avg, "1.0" );	        /* Averaging */
    
/* the main menu loop */

    while ( true )
    {
        item = ReadMenuItem();
        switch ( item )
        {
            case item_EXEC:
		DoExec();
                break;
            case item_SET:
		ReadOption();
		batch = OptionTrue( opt_batch );
                if ( !Check_General_Options () ) {
	            printf( "Error in main: general options set are not okay\n" );
	            MError( );
                }
                terseness   = GetOptionIVal ( opt_terseness );                
                break;
            case item_range:
		oldmenu = SetOneMenu( "Give name of ranges file (or - to end):" );
		while ( !done )
		{
		    MReadWord( rng_filnam );
		    if ( *rng_filnam == '-' || strlen( rng_filnam ) == 0 )
			break;
		    if ( ! ( done = OPENr( fp, rng_filnam ) ) )
		    {
			printf( " * ranges file <%s> could not be opened\n", rng_filnam );
			MError();
		    }
		}
		ResetMenu( oldmenu );
		if ( done ) 
		{
		    fclose( fp );
		    range = ReadRanges( rng_filnam );
		}
                break;		
            case item_do:
		goto doit;
            case item_QUIT:
            case item_NULL:
                exit( 0 );
                break;
            default:
                ProcessMenuError( menu_fatal );
                break;
        }
    }

/* do the main thing */

doit:

/* process the distance-values file saving the header lines and file names */

    *dval_filnam = '\0';
    ProcessDistvalFile( &val_list, &dist_count, &struc_count, &models, &qtv, 
                                    theSet.filnam, dinf_filnam, drst_filnam,
		                    ident, true, &ftype, dval_filnam );
    if ( ftype != inf_extended )
    {
	FATAL( " * Wrong file/data type - Code error - case 4" );
    }
/*
    else
    {
	dist_val  = NULL;
	dist_vall = val_list.lst_desc->data;
    }
*/

/* process the distance-restraints file (file name read from distance-values file) */

    if ( IsEmptyStr( drst_filnam ) )
    {
     FATAL( " * Name of dist-restr file not found on distance-values file");
    }
    dist_rst = ProcessDistrstFile( &count, &qtype, filnam_s, filnam_i, 
                                   ident, false, drst_filnam );
    if ( NOTEQ( theSet.filnam, filnam_s ) ||
         NOTEQ( dinf_filnam, filnam_i )   ||
         dist_count != count )
    {
        FATAL(" * Inconsistent information on values and restraints files");
    }

/* process the distance-info file (file name read from distance-values file) */


    ProcessDistinfFile( &inf_list, &count, &s_count, &qtype, filnam_s, filnam_r,
                                   ident, false, &ftype, dinf_filnam );
    if ( ftype != inf_extended )
    {
	FATAL( " * Wrong file/data type - Code error -case 5" );
    } else
    {
	/*dist_inf  = inf_list.lst_desc->data;*/
	dist_infl = (Distinfl *) (inf_list.lst_desc->related->data);
    }

    if ( NOTEQ( theSet.filnam, filnam_s ) ||
         NOTEQ( dinf_filnam, filnam_i )   ||
         dist_count != count )
    {
        FATAL( " * Inconsistent information on values and info files" );
    }

    if ( terseness < 6 ) {
        printf ("Read %d restraints\n", dist_count);
    }


/* open and read the strucset file */
    if ( terseness < 6 ) {
        printf("Open and read the strucset file\n");
    }

    if ( EMPTYSET( theSet ) )
    {
	if ( !ProcessStrucsetFile( &theSet, -1 ) )
	    exit( 1 );
    }     
    
    molptr = theSet.molecule;
    
/* AquaCompleteness specific routines */


/* Processing the options */

    save_pseudos    = OptionTrue( opt_savpseu );
    del_red_l       = OptionTrue( opt_delredl );
    del_red_u       = OptionTrue( opt_delredu );
    use_distance    = OptionTrue( opt_usedist );

    atominp         = GetOptionCVal( opt_atominp );
    thresh          = GetOptionRVal( opt_thresh );
    avg             = GetOptionRVal( opt_avg );


/* Check the limits of the options */  

    if ( !Check_Options () ) {
	nrerror( "Error in main: options are not okay");
    }

    if ( strlen( rng_filnam ) == 0 ) {
	printf("opt_atominp is empty\n");	     
	nrerror( "Error in main: options are not okay");
    }
    
/* Printing the options */
    if ( terseness < 6 ) {

        printf ("\nThe following OPTIONS are used\n");

        printf ("%8d shells for pairs\n",		GetOptionIVal( opt_numb_shells_prs ) );
        printf ("%8.3f min distance pairs\n",	GetOptionRVal( opt_min_dist_prs )    );
        printf ("%8.3f max distance pairs\n\n",	GetOptionRVal( opt_max_dist_prs )    );

        printf ("%8d shells for restraints\n",      GetOptionIVal( opt_numb_shells_rst ) );
        printf ("%8.3f min distance restraints\n",	GetOptionRVal( opt_min_dist_rst )    );
        printf ("%8.3f max distance restraints\n\n",GetOptionRVal( opt_max_dist_rst )    );

        printf ("Observable atoms from set                    : %s\n", 
		    atominp);
        printf ("Use intra-residual contacts                  : %c\n", 
		    Boolean_To_Char( OptionTrue( opt_use_intra ) ) );
        printf ("Threshold redundancy                         : %-9.2f\n", 
		    thresh );
        printf ("Do not check pseudos for redunancy           : %c\n", 
		    Boolean_To_Char( OptionTrue( opt_savpseu ) ) );
        printf ("Count restraints by atom                     : %c\n", 
		    Boolean_To_Char( OptionTrue( opt_double_count ) ) );
        printf ("Show some extra detailed information         : %c\n", 
		    Boolean_To_Char( OptionTrue( opt_detail ) ) );
        printf ("Averaging power                              : %-9.2f\n", 
		    avg );
        printf ("Range file (- for none)                      : %s\n", 
		    rng_filnam );

        printf ("Analyse the SSA                              : %c\n",
	        Boolean_To_Char( OptionTrue( opt_anal_ssa ) ) );
    }
            
/* MATCH TRIPLETS PRESENT IN THE MOL WITH SSA RESTRAINTS
    WRITE RESULTS TO TEMPORARY FILE  */
    
    if ( OptionTrue( opt_anal_ssa ) ) {
        if ( terseness < 6 ) {
	    printf("Analyse the theoretical and experimental stereospecific info\n");
        }
	Match_SSA ( &theSet, dist_infl, dist_count );
    }

    if ( terseness < 6 ) {
        printf("Open and read the strucset file\n");
    }
    
    sprintf(theSet_ref.filnam, "%s/%s.set", PathName, atominp);
   
    if ( EMPTYSET( theSet_ref ) )
    {
	if ( !ProcessStrucsetFile( &theSet_ref, -1 ) )
	    exit( 1 );
    }
    

/* MAKES OBSERVABLE ATOMLIST (ONLY FIRST MOL) */

    obs_atom_count = Obs_Atoms ( &obs_atom, &theSet, &theSet_ref );     

    if ( obs_atom == NULL ) {
        printf("No observable atoms !!\n");
	nrerror( "Error in main: Wrong init observable atoms" );
    }



/* MAKE "OBSERVABLE" RESTRAINT LIST */
    
    if ( terseness < 6 ) {
        printf ("Restraints can be collapsed to become \"observable\"\n");
    }
    /* expansie mogelijk maken door dist_infl uit te breiden met qhow te maken */


    extra_distinfl = CREATE_NEW( Extra_Distinfl, dist_count );

    /* Only fill the cor_A and cor_B values !!! */
    Calculate_Corrections ( extra_distinfl, dist_infl, dist_count);
    
    obs_restr_count = Obs_Restraints ( 
		    &obs_restr, dist_infl, dist_rst,dist_count, obs_atom, &theSet, 
		    &rev_mapping_table, extra_distinfl, 
		    &number_not_observable, &number_not_in_range, &number_ssa_double,     
		    &number_intra,          &number_low_only,     &number_nobs_atom );
    
    if ( obs_restr == NULL ) {
	nrerror( " * Wrong init observable restraints" );
    }
        
    
    DR_list.list = obs_restr;
    DR_list.dist_count = obs_restr_count;
    DR_list.struc_count = 1;

/* multiples uit restraint lijst verwijderen ??? */
    
    multi_restr_count = CleanupDRList( &DR_list, 0, &exit_flag );
    
    if ( terseness < 9 ) {
        printf ("Found %d multiply occuring restraints.\n", multi_restr_count);
        printf ("Leaving %d restraints.\n", obs_restr_count  
                                            - multi_restr_count);
    }


/* Mark redundant (only fixed ones) restraints with keep flag false 
    using the AquaRedun subroutine: */

    if ( OptionTrue( opt_use_intra ) ) {
	red_restr_count = Check_Redundant (obs_restr, obs_restr_count, multi_restr_count, 
			     OptionTrue( opt_detail ), 
			     save_pseudos, del_red_l, del_red_u, use_distance, thresh );
        if ( terseness < 9 ) {
	    printf ("Marked %d redundant restraints as such.\n", red_restr_count);
        }
    } else {
	red_restr_count = 0;
    }
    
    no_restr_count = multi_restr_count + red_restr_count;
    if ( terseness < 9 ) {
        printf("Leaving %d restraints.\n", obs_restr_count - no_restr_count );
        printf("Writing the set of %d observable restraints.\n", 
	obs_restr_count-no_restr_count);
    }

    ProcessDROutput( &DR_list, no_restr_count, 
		     "noe", false, molptr, struc_count);

/* MAKE A LIST OF PAIRS and a list of distances for each in parallel
 * that have an average distance smaller than the CUTOFF.
 */


    obs_pair_count = Obs_Pair_Val ( &pair, &pair_val, obs_atom, obs_atom_count, &theSet);

				      
    if ( pair == NULL ) {
	nrerror( " * Wrong init observable pairs" );
    }

    DR_list.list = pair;
    DR_list.dist_count = obs_pair_count;
    DR_list.struc_count = 1;

/* Mark redundant (only fixed ones) pairs with keep flag false 
    using the AquaRedun subroutine: */
    

    if ( OptionTrue( opt_use_intra ) ) {
    
	red_pair_count = Check_Redundant  ( pair, obs_pair_count, 0, 
			    OptionTrue( opt_detail ), 
			    save_pseudos, del_red_l, del_red_u, 
			    use_distance, 
			    (float) 9999.9 );
	
        if ( terseness < 9 ) {
	    printf ("Marked %d redundant pairs as such.\n", red_pair_count);
        }
    } else {
	red_pair_count = 0;
    }
    no_pair_count = red_pair_count;
    
    if ( terseness < 9 ) {
        printf ("Leaving %d restraints.\n", obs_pair_count - no_pair_count);
    }

    if ( terseness < 9 ) {
        printf("Writing the set of %d observable pairs.\n", 
	    obs_pair_count - no_pair_count);
    }
    ProcessDROutput( &DR_list, no_pair_count,"noe",false,molptr,struc_count);
        
    overall_completeness_factor = Process_Distances_Shells ( 
	atominp, dist_count, theSet, theSet_ref, obs_atom, obs_atom_count,
	obs_restr, obs_restr_count, no_restr_count, 
	pair, pair_val, obs_pair_count,  no_pair_count, 
	drst_filnam, filnam, dist_infl, rev_mapping_table, 
	rng_filnam );
 
    
    printf("OVERALL COMPLETENESS FACTOR is %8.3f or %d percent\n", 
     overall_completeness_factor, (int) (100*overall_completeness_factor + 0.5) );

    if ( OptionTrue( opt_anal_ssa ) ) {
	if ( !Concatenate_file_A_To_B ( "AquaCompl.temp", filnam ) ) {
	    nrerror("Couldn't merge the two result files.");
	}	
    }

    if ( terseness < 6 ) {
        printf("Finished %s, have a nice day\n\n", program_name);
    }
    return(0);
}





static int Obs_Atoms ( Boolean **obs_atom, 
			Strucset *setptr, Strucset *setptr_ref )          

/* Return a Boolean list of all atoms in theSet.molecule F/T observable */

{

    int	    i, j, count, idx;
    Boolean selected;
    Residue *residues;
    Atom    *atoms, *atoms_ref;
    Mol	    *molptr, *molptr_ref;
    
    molptr	= setptr->molecule;
    molptr_ref  = setptr_ref->molecule;

    atoms	= molptr->atoms;
    atoms_ref	= molptr_ref->atoms;

    residues	= molptr->residues;

    *obs_atom = CREATE_NEW( Boolean, molptr->atom_count);
    
    count = 0;
    
    for (i=0; i<molptr->atom_count; i++) {
    
	(*obs_atom)[i] = false;
	
	if ( range ) {
	    idx = FindResidueIndexFromNum( atoms[i].res_num_new, molptr );
	    selected = idx < 0 ? false : InResidueRange( residues+idx );
	    if ( !selected ) {
		continue;
	    }
	}
	
	for (j=0; j<molptr_ref->atom_count; j++) {
	    
	    if ( EQID ( atoms[i].atom_int, atoms_ref[j].atom_int ) &&
		 EQID ( atoms[i].res_nam,  atoms_ref[j].res_nam ) ) {
		
		(*obs_atom)[i] = true;
		count++;
		break;
	    }
	}
    }
    
    if ( ( count < 0 ) || ( count > molptr->atom_count ) ) {
	nrerror( "Error in Obs_Atoms: bad count in function Obs_Atoms");
    }
    
    if ( terseness < 9 ) {
        printf ("Found %5d total number atoms in structure set\n", molptr->atom_count);
        printf ("Found %5d \"observable\" atoms in structure set\n", count);
    }
    
    return ( count );
}



static int Obs_Restraints ( DRestraint **DR, 
    Distinfl *dist_infl, Distrst *dist_rst, int dist_count,
    Boolean *obs_atom, Strucset *setptr, 
    int **rev_mapping_table, Extra_Distinfl *extra_distinfl,
    int *number_not_observable, int *number_not_in_range, int *number_ssa_double,     
    int *number_intra,          int *number_low_only,     int *number_nobs_atom ) 

{

/*  Make a list of restraints from the (collapsed) atoms of the
    original restraints */    

    int		i, j, count;
    int         j_found = (int) NOT_A_NUMBER ;
    int		number_expandable, number_collapsed;
    int		idx_A, idx_B;
    int		idx_P_A, idx_P_B, idx_Q_A, idx_Q_B;
    Boolean	found_a, found_b, expand_a, expand_b, collapse_a, collapse_b;
    Boolean	*restraint_in_range, *ssa_double, select_A, select_B;
    Boolean	use_intra, detail, found;
    Distval     correction;
    FILE	*fp=stdout;
    Mol		*molptr;
    Atom	*atoms;
    Residue	*residues;
    
    number_collapsed        = 0;
    number_expandable       = 0;
    
    *number_not_observable  = 0;
    
    *number_not_in_range    = 0;
    *number_ssa_double      = 0;
    *number_low_only        = 0;
    *number_intra           = 0;
    *number_nobs_atom	    = 0;
    
    count = 0;
    
    molptr	= setptr->molecule;
    atoms	= molptr->atoms;
    residues	= molptr->residues;

    use_intra = OptionTrue( opt_use_intra );
    detail    = OptionTrue( opt_detail );
    
    *DR = CREATE_NEW( DRestraint, dist_count ); 
    *rev_mapping_table = CREATE_NEW( int, dist_count );
    Fill_Int_Array ( (*rev_mapping_table), -1, 0, dist_count -1 ); 
    restraint_in_range = CREATE_NEW( Boolean, dist_count );
    ssa_double = CREATE_NEW( Boolean, dist_count );
    
    if ( terseness < 6 ) {
        printf("FIRST deleting the restraints that are not in the range\n");
    }
        
    if ( range ) {
	for (i=0; i<dist_count; i++) {
	    idx_A = FindResidueIndexFromNum( 
		atoms[ dist_infl[i].A[0].idx ].res_num_new, molptr );
	    idx_B = FindResidueIndexFromNum( 
		atoms[ dist_infl[i].B[0].idx ].res_num_new, molptr );
	    select_A = idx_A < 0 ? false : InResidueRange( residues+idx_A );
	    select_B = idx_B < 0 ? false : InResidueRange( residues+idx_B );
	    if (  select_A && select_B ) {
		restraint_in_range[i] = true;
	    } else {
		restraint_in_range[i] = false;
		(*number_not_in_range)++;
		if  ( detail ) {
		    printf("Not both in range (%4d) ", i);
		    idx_A = dist_infl[i].ass_A;
		    idx_B = dist_infl[i].ass_A;
		    WriteFAtom( fp, &(atoms[ dist_infl[i].A[ idx_A ].idx ]), use_int );
		    WriteFAtom( fp, &(atoms[ dist_infl[i].B[ idx_B ].idx ]), use_int );
		    printf("\n");
		}
	    }
	}

        if ( terseness < 9 ) {
	    printf("Found %d restraints not both in range; so not observable\n", 
	    *number_not_in_range );	
        }

    } else {
	Fill_Boolean_Array( restraint_in_range, true, 0, dist_count-1 );
    }

    
    if ( terseness < 9 ) {
        printf("SECOND Deleting the restraints with a pseudoatom for which \n");
        printf("the same restraints with a stereospecific assigned atom exists\n");
    }
        
    /* Procedure will also delete normal double restraints, if present */
		
    Fill_Boolean_Array( ssa_double, false, 0, dist_count-1 );		
		
    /* Loop over the to check restraints P with counter i */		
    for (i=0; i<dist_count; i++) {
	
	if ( ! restraint_in_range[i] ) {
	    continue;
	}

	found = false;

	/* FIRST ATOM  */

	if ( dist_infl[i].ass_A == 2 ) {
	    idx_P_A = dist_infl[i].A[ 0 ].idx;

	    if ( dist_infl[i].ass_B == 2 ) {
		idx_P_B = dist_infl[i].B[ 0 ].idx;
	    } else {
		idx_P_B = dist_infl[i].B[ dist_infl[i].ass_B ].idx;
	    }
	    
	    for (j=0; j<dist_count; j++) {
    
		if ( ( i==j ) || ! restraint_in_range[j] ){
		    continue;
		}
    
		idx_Q_A = dist_infl[j].A[ 0 ].idx;
		
		if ( dist_infl[i].ass_B == 2 ) {
		    idx_Q_B = dist_infl[j].B[ 0 ].idx;
		} else {
		    idx_Q_B = dist_infl[j].B[ dist_infl[j].ass_B ].idx;
		}
    
		if ( ( ( idx_P_A == idx_Q_A ) && ( idx_P_B == idx_Q_B ) ) ||
		     ( ( idx_P_A == idx_Q_B ) && ( idx_P_B == idx_Q_A ) ) ) {
		    found = true;
		    j_found = j;
		    break;
		}
	    }
	}
	
	/* SECOND ATOM */

	if ( ( !found ) && ( dist_infl[i].ass_B == 2 ) ) {
	    
	    idx_P_B = dist_infl[i].B[ 0 ].idx;

	    if ( dist_infl[i].ass_A == 2 ) {
		idx_P_A = dist_infl[i].A[ 0 ].idx;
	    } else {
		idx_P_A = dist_infl[i].A[ dist_infl[i].ass_A ].idx;
	    }
	    
	    for (j=0; j<dist_count; j++) {
    
		if ( ( i==j ) || ! restraint_in_range[j] ){
		    continue;
		}
    
		idx_Q_B = dist_infl[j].B[ 0 ].idx;
		
		if ( dist_infl[i].ass_A == 2 ) {
		    idx_Q_A = dist_infl[j].A[ 0 ].idx;
		} else {
		    idx_Q_A = dist_infl[j].A[ dist_infl[j].ass_A ].idx;
		}
    
		if ( ( ( idx_P_A == idx_Q_A ) && ( idx_P_B == idx_Q_B ) ) ||
		     ( ( idx_P_A == idx_Q_B ) && ( idx_P_B == idx_Q_A ) ) ) {
		    found = true;
		    j_found = j;
		    break;
		}
	    }

	}
		
	if ( found ) {
    
	    ssa_double[i] = true;
	    (*number_ssa_double)++;
		    
	    if  ( detail && ( terseness < 9 )) {
		j =  j_found;

		printf("SSA double restraint (%4d with %4d) ", i, j);

		printf("[%4d %4d ** %4d %4d] ", 
		    dist_infl[i].A[ dist_infl[i].ass_A ].idx, 
		    dist_infl[i].B[ dist_infl[i].ass_B ].idx, 
		    dist_infl[j].A[ dist_infl[j].ass_A ].idx, 
		    dist_infl[j].B[ dist_infl[j].ass_B ].idx);

		WriteFAtom( fp, &(atoms[ dist_infl[i].A[ dist_infl[i].ass_A ].idx ]), use_int );
		WriteFAtom( fp, &(atoms[ dist_infl[i].B[ dist_infl[i].ass_B ].idx ]), use_int );

		printf(" ****** ");

		WriteFAtom( fp, &(atoms[ dist_infl[j].A[ dist_infl[j].ass_A ].idx ]), use_int );
		WriteFAtom( fp, &(atoms[ dist_infl[j].B[ dist_infl[j].ass_B ].idx ]), use_int );

		printf("\n");
	    }
	}
    }

    if ( terseness < 9 ) {
        printf ("Found %d restraints SSA double; so not observable\n", 
	    *number_ssa_double );	
    }


    for (i=0; i<dist_count; i++) {		

	if ( !restraint_in_range[i] || ssa_double[i] ) {
	    continue;
	}

	if ( dist_rst[ i ].typ_bnd == 'l') {
	    if  ( detail && ( terseness < 9 ) ) {
		printf("lowerbound only,  is not usable for restraint %d\n", i);
	    }
	    (*number_low_only)++;
	    continue;
	}

	if ( ! use_intra ) {
	    if ( Same_Residue ( atoms[ dist_infl[i].A[0].idx ],  
	                        atoms[ dist_infl[i].B[0].idx ] ) ) {
		if  ( detail && ( terseness < 9 ) ) {
		    printf("Excluded intraresidual (%4d) ", i);
		    WriteFAtom( fp, &(atoms[ dist_infl[i].A[0].idx ]), use_int );
		    WriteFAtom( fp, &(atoms[ dist_infl[i].B[0].idx ]), use_int );
		    printf("\n");
		}
		(*number_intra)++;
		continue;
	    }
	}


	found_a = false;
	found_b = false;
	expand_a =  false;
	expand_b =  false;
	collapse_a =  false;
	collapse_b =  false;


	/* Checking partner A */

	if ( obs_atom[  dist_infl[i].A[ dist_infl[i].ass_A ].idx  ] ) {
	    found_a = true;
	}
	
	if ( !found_a ) {               /* Can it be expanded/collapsed ? */
	
	    if ( dist_infl[i].count_A != 1 ) {

					/* Does it then become observable ?*/
    
		if ( dist_infl[i].ass_A == 2 ) {              /* Expandable */
    
		    if ( obs_atom[ dist_infl[i].A[0].idx ] ||
		         obs_atom[ dist_infl[i].A[1].idx ]    ) {
			found_a = true;
			expand_a = true;
			if ( detail && ( terseness < 6 ) ) {
			    printf("Expandable partner a of restraint %d\n", i);
			}
		    }
    
		} else {                                    /* Collapse it */
    
		    if ( obs_atom[ dist_infl[i].A[2].idx ] ) {
			found_a = true;
			collapse_a = true;
			if ( detail && ( terseness < 6 )) {
			    printf("Collapsing partner a of restraint %d\n", i);
			}
		    }
		}
	    }
	}	
	
	

	/* Checking partner B */

	if ( obs_atom[  dist_infl[i].B[ dist_infl[i].ass_B ].idx  ] ) {
	    found_b = true;
	}
	
	if ( !found_b ) {               /* Can it be expanded/collapsed ? */
	
	    if ( dist_infl[i].count_B != 1 ) {

					/* Does it then become observable ?*/
    
		if ( dist_infl[i].ass_B == 2 ) {              /* Expandable */
    
		    if ( obs_atom[ dist_infl[i].B[0].idx ] ||
		         obs_atom[ dist_infl[i].B[0].idx ]    ) {
			found_b = true;
			expand_b = true;
			if ( detail && ( terseness < 6 )) {
			    printf("Expandable partner b of restraint %d\n", i);
			}
		    }
    
		} else {                                    /* Collapse it */
    
		    if ( obs_atom[ dist_infl[i].B[2].idx ] ) {
			found_b = true;
			collapse_b = true;
			if ( detail && ( terseness < 6 )) {
			    printf("Collapsing partner b of restraint %d\n", i);
			}
		    }
    
		}
	    }
	}	
	
	if ( collapse_a || collapse_b || expand_a || expand_b) {
	    /*
		printf("Collapsed or expandable: ");
		WriteFAtom( fp, &(atoms[ dist_infl[i].A[ dist_infl[i].ass_A ].idx ]), use_int );
		WriteFAtom( fp, &(atoms[ dist_infl[i].B[ dist_infl[i].ass_B ].idx ]), use_int );
		printf("\n");
	    */
	}
	
	if ( found_a && found_b ) {
	    
	    /* Collapsed restraint (from stereospecific to pseudo) can be
	     * generated more than once so the list should later on be 
	     * corrected for doubles.
	     */
	    	    
	
	    /* Just count the expandable restraints DON't actually expand */
	    
	    if ( expand_a || expand_b ) {
		number_expandable++;
	    }
	    
	    if ( collapse_a || collapse_b ) {
		number_collapsed++;
	    }
	    
	                        /* Add the (collapsed) restraint(s) */

	    correction = 0.0;
	    
	    if ( collapse_a ) {
		idx_A = 2;
		if ( Same_Floats(   (Distval) (extra_distinfl[i].cor_A), 
				    (Distval) BAD_CORRECT ) ) {
		    nrerror("Error in Obs_Restraints A");
		}
		correction += (Distval) extra_distinfl[i].cor_A;
	    } else {
		idx_A = dist_infl[i].ass_A;
	    }			

	    if ( collapse_b ) {
		idx_B = 2;
		if ( Same_Floats(   (Distval) (extra_distinfl[i].cor_B), 
				    (Distval) BAD_CORRECT ) ) {
		    nrerror("Error in Obs_Restraints B");
		}
		correction += (Distval) extra_distinfl[i].cor_B;
	    } else {
		idx_B = dist_infl[i].ass_B;
	    }			
	
	    (*DR)[count].idx_A = dist_infl[i].A[ idx_A ].idx ;
	    (*DR)[count].atom_A = *(atoms + (*DR)[count].idx_A );
    
	    (*DR)[count].idx_B = dist_infl[i].B[ idx_B ].idx;
	    (*DR)[count].atom_B = *(atoms + (*DR)[count].idx_B );
	    	    
	    (*DR)[count].keep   = true;
	    
	    switch ( dist_rst[i].typ_bnd ) {
	
		case 'u' :  (*DR)[count].stat   = bnd_upper;
			    (*DR)[count].lo_bnd = 0.0;
			    (*DR)[count].up_bnd = dist_rst[i].up_bnd + correction;
			    (*DR)[count].dist   = 0.0;
			    break;
		case 'l' :  nrerror("ERROR in Obs_Restraints: wrong case");
			    break;
		case 'b' :  (*DR)[count].stat   = bnd_uplo;
			    (*DR)[count].lo_bnd = dist_rst[i].lo_bnd - correction;
			    (*DR)[count].up_bnd = dist_rst[i].up_bnd + correction;
			    (*DR)[count].dist   = 0.0;
			    break;
		case 'B' :  (*DR)[count].stat   = bnd_uld;
			    (*DR)[count].lo_bnd = dist_rst[i].lo_bnd - correction;
			    (*DR)[count].up_bnd = dist_rst[i].up_bnd + correction;
			    (*DR)[count].dist   = dist_rst[i].dist;
			    break;
		default  :  nrerror("ERROR in Obs_Restraints: default case");
	    }
	    
	    (*rev_mapping_table)[count++] = i;
	    	    
/*
	    if ( DEBUG_INFO ) {
		printf("%d %d %d %d %8.3f\n", count, i, idx_A, idx_B, correction); 
	    }
*/
	} else {
	
	    (*number_nobs_atom)++;
	    
	    printf("Restraint (%4d) not observable ", i);
	    WriteFAtom( fp, &(atoms[ dist_infl[i].A[0].idx ]), use_int );
	    WriteFAtom( fp, &(atoms[ dist_infl[i].B[0].idx ]), use_int );
	    printf("\n");
	}
	
    }
    
    /*
    for (i=0; i<dist_count; i++) {
	printf("Map %d %d\n", i, (*rev_mapping_table)[i]);
	if ( (*rev_mapping_table)[i] < -1 ) {
	    nrerror("Error in Obs_Restraints");
	}
    }
    */
    
    *number_not_observable += 
	*number_not_in_range + *number_ssa_double +
        *number_low_only     + *number_intra +      *number_nobs_atom;
			    
    if ( dist_count != *number_not_observable + count ) {
	printf (" %d + %d != %d restraints \n", *number_not_observable, count, dist_count);
	nrerror("Error in Obs_Restraints; bad counting the restraints");
    }

    printf ("\n\nFound %d restraints in file.\n", dist_count);
    printf ("Found %5d restraints not observable.\n", *number_not_observable);
    printf ("-> Excluded %5d restraints not in range\n",             *number_not_in_range);
    printf ("   Excluded %5d ssa double restraints\n",		     *number_ssa_double);
    printf ("   Excluded %5d low only restraints\n",		     *number_low_only);
    printf ("   Excluded %5d intraresidual restraints\n",            *number_intra);
    printf ("   Excluded %5d restraints with not observable atoms\n",*number_nobs_atom);
    printf ("Now   %5d restraints \"observable\".\n", count);

    printf ("Collapsed %d restraints\n", number_collapsed );
    printf ("Found %d expandable restraints\n\n\n", number_expandable );
    
    return ( count );
}

    
static int Obs_Pair_Val ( DRestraint **pair, Distvall **distances, 
                      Boolean *obs_atom, int obs_atom_count,
		      Strucset *setptr )
{

/* Returns a list of pairs and values that have an average distance below
 * the cutoff.
 */

#define MAX_PAIRS 40000

    int obs_pair_count;                    /* Will be returned */
    int i, j, k /*l*/ ;
    int max_obs_pairs, obs_atom_count_act;
    /*FILE *fp=stdout;*/
    Boolean use_intra;
    Distval *distval, lo, up, average;
    Distval avg, avg_div;
    Distval cut_distance, sum;
    Atom *atoms;
    Mol *molptr;
    Struc *strucptr;
    
    molptr = setptr->molecule;
    atoms  = molptr->atoms;
    
    obs_pair_count = 0;
    obs_atom_count_act = 0;

    use_intra = OptionTrue( opt_use_intra );
    cut_distance = GetOptionRVal( opt_max_dist_prs );
    avg = GetOptionRVal( opt_avg );
    avg_div = (Distval) 1.0 / avg;
    
    distval = CREATE_NEW( Distval, setptr->struc_count ); 
           
    
    /* The values (most expensive) themselve are located dynamically */
    max_obs_pairs = (int) ( 0.5 * obs_atom_count * (obs_atom_count - 1) + 0.5 );

    if ( max_obs_pairs > MAX_PAIRS ) {
	printf("WARNING in Obs_Pair_Val: too many pairs\n ");
	printf("MAY LEAD TO ABORT WHEN TOO MANY ARE FOUND\n");
	max_obs_pairs = MAX_PAIRS;
	printf("maximum number of pairs is set to %d\n", MAX_PAIRS);
    }
    
    *pair = CREATE_NEW( DRestraint, max_obs_pairs ); 

    *distances = CREATE_NEW( Distvall, max_obs_pairs ); 
           
    for (i=0; i<molptr->atom_count; i++) {
	
	if ( !obs_atom[i] ) {
	    continue;
	}
	
	obs_atom_count_act++;
	
	for (j=i+1; j<molptr->atom_count; j++) {
	    
	    if ( !obs_atom[j] ) {
		continue;
	    }
	    
	    if ( ! use_intra ) {
		if ( Same_Residue (atoms[i], atoms[j]) ) {
		    continue;
		}
	    }
	    
	    sum = 0.0;
	    strucptr = setptr->fst_struc;
	    
	    for (k=0; k<setptr->struc_count; k++) {

		distval[k] = sqrt( SQRDIST( strucptr->coords[ i ], 
		                            strucptr->coords[ j ] ) );
								
		sum  += pow( distval[k], avg );
		strucptr = strucptr->next_struc;
	    }
	    
	    average = pow ( sum / setptr->struc_count, avg_div );
	    
	    if ( average > cut_distance ) {
		continue;	    
	    }
	    
	    
	    /*
	    if ( DEBUG_INFO ) {
		WriteFAtom( fp, &(setptr->molecule->atoms)[ i ], use_int );
		WriteFAtom( fp, &(setptr->molecule->atoms)[ j ], use_int );
		for (l=0; l<k; l++) {
		    printf("%9.3f ", distval[l]);
		}
		printf("\n");
	    }
	    */
	    
	    if (obs_pair_count > (max_obs_pairs-1) ) {
		nrerror("Error in Obs_Pair_Val: found too many pairs");
	    }
	    
	    /* Find out the lo, up and average of the array */
	    
	    lo =      Lowest_Array  ( distval, setptr->struc_count );
	    up =      Highest_Array ( distval, setptr->struc_count );
	    
	    (*pair)[obs_pair_count].idx_A  = i;
	    (*pair)[obs_pair_count].atom_A = *(atoms + i );
    
	    (*pair)[obs_pair_count].idx_B  = j;
	    (*pair)[obs_pair_count].atom_B = *(atoms + j );
    
	    (*pair)[obs_pair_count].keep   = true;
	    (*pair)[obs_pair_count].stat   = bnd_uld;

	    (*pair)[obs_pair_count].lo_bnd = lo;
	    (*pair)[obs_pair_count].up_bnd = up;
	    (*pair)[obs_pair_count].dist   = average;
	    
	    
		    /* Add the distances (dynamically) for this pair */
				    
	    (*distances)[ obs_pair_count ].values = 
				    CREATE_NEW( Distval, setptr->struc_count ); 

	    (*distances)[ obs_pair_count ].count = setptr->struc_count;

	    Copy_Array ( distval, (*distances)[obs_pair_count].values, 
			setptr->struc_count );
/*	    
	    if ( DEBUG_INFO ) {
		for (k=0; k<setptr->struc_count; k++) {
		    printf ("%8.3f ", (*distances)[obs_pair_count].values[k]);
		}
		printf ("\n");
	    }
*/
	    obs_pair_count++;					    		    		    
	}
    }
    

    if ( obs_atom_count != obs_atom_count_act ) {
	printf ("There are two different counts of \"observable\" atoms: %d %d\n", 
	          obs_atom_count, obs_atom_count_act);
	nrerror("Error in Obs_Pair");
    }

    if ( terseness < 6 ) {
        printf ("Found %5d  \"observable\" pairs.\n", obs_pair_count);
    }
    
    return ( obs_pair_count );
}




static Distval Lowest_Array ( Distval *values, int count ) 
{    
    int i;
    Distval dist;
    
    i = 1;
    dist = values[0];  

    if ( count == 0 ) {
	nrerror("Error in Lowest_Array: count is zero");
    }

    while ( i < count ) {
	if ( values[i] < dist ) {
	    dist = values[i];
	}
	i++;
    }
    return ( dist );
}



static Distval Highest_Array ( Distval *values, int count ) 
{    
    int i;
    Distval dist;
    
    i = 1;
    dist = values[0];  
    
    if ( count == 0 ) {
	nrerror("Error in Highest_Array: count is zero");
    }

    while ( i < count ) {
	if ( values[i] > dist ) {
	    dist = values[i];
	}
	i++;
    }
    return ( dist );
}



static Distval Power_Average_Distval_Array ( Distvall distvall, Distval avg ) 
{    
    int i;
    Distval dist, avg_div, sum;
    
    i = 0;
    sum = 0.0;
    avg_div = 1.0 / avg;

    if ( distvall.count == 0 ) {
	nrerror("Error in Power_Average_Distval_Array: count is zero");
    }
        
    while ( i < distvall.count ) {
	sum += pow( distvall.values[i], avg );
	i++;
    }
    dist = pow ( sum / distvall.count, avg_div );
    
    return ( dist );    
}



static void Copy_Array ( Distval *a, Distval *b, int count )
{
    /* Copy array a to b */
    
    int i;
    
    for (i=0;i<count;i++) {
	b[i] = a[i];
    }
    
    return;
    
}
    

static Distval Process_Distances_Shells ( 
    char *atominp, int dist_count, 
    Strucset theSet, Strucset theSet_ref, Boolean *obs_atom, int obs_atom_count,
    DRestraint *obs_restr,                int obs_restr_count, int no_restr_count, 
    DRestraint *pair, Distvall *pair_val, int obs_pair_count,  int no_pair_count, 
    Filnam drst_filnam, Filnam filnam, 
    Distinfl *dist_infl, int *rev_mapping_table,  
    Filnam rng_filnam )
{

    Distval total_compl_fac, compl_fac, maxdist_prs, bin_prs;
    int	    *S_pair_ass;
    Shells  *S_T_restr_ass;
    int	    j, s;
    int	    *pair_shell_count, shell_count;
    int	    pair_total_count, shell_total_count;
    Boolean *pair_matched;

    DRList      DR_list;
    
    pair_shell_count = CREATE_NEW( int, GetOptionIVal( opt_numb_shells_prs )+1);

    pair_matched = CREATE_NEW(Boolean, obs_pair_count);
    for (j=0; j< obs_pair_count; j++) {
	pair_matched[j] = false;
    }
    
    if ( terseness < 9 ) {
        printf("\nProcessing the calculated distances \n");
        printf("and match them to the observed restraints\n\n\n");
    }
    
    /* Storing the info on classification in two sets of (doubly) int's 
	one for the pairs and one for the restraints
	watch out that the list run in parallel with redundant/multiple 
	restraints/pairs marked with a false keep flag
    */

    shell_total_count = 0;
    pair_total_count  = 0;

    S_pair_ass = CREATE_NEW( int, obs_pair_count );

    Fill_Int_Array( S_pair_ass,  0, 0, obs_pair_count-1 );
    
    S_T_restr_ass = CREATE_NEW( Shells, obs_restr_count );

    Fill_2Int_Array( S_T_restr_ass,  0, 0, 0, obs_restr_count-1 );
    
    
/* Loop over expanding shell S */
    
    bin_prs = ( GetOptionRVal( opt_max_dist_prs ) - 
		GetOptionRVal( opt_min_dist_prs ) ) / 
		GetOptionIVal( opt_numb_shells_prs );

    maxdist_prs =    GetOptionRVal( opt_min_dist_prs );
    
    for (s=1; s<=GetOptionIVal( opt_numb_shells_prs ); s++) {

	maxdist_prs += bin_prs;	    

	/* Categorize pairs for this shell, set S_pair */
	
	pair_shell_count[s] = Categorize_Pairs ( maxdist_prs, s, 
	    pair, pair_val, S_pair_ass, obs_pair_count, no_pair_count );
	
	pair_total_count += pair_shell_count[s];
		
        if ( terseness < 6 ) {
	    printf("%8d pairs in shell %d until %8.3f Angstrom\n", 
	        pair_shell_count[s], s, maxdist_prs);
        }
    }
    
    if ( pair_total_count != obs_pair_count - no_pair_count ) {
	printf("Found %8d pair_total_count\n", pair_total_count);
	nrerror("Error in Process_Distances_Shells: pair count problem");
    }    
    
/* Again loop over expanding shell S */

   for (s=1; s<=GetOptionIVal( opt_numb_shells_prs ); s++) {
   
	/* Match observable distances against restraints 
	    and check in which t bin they fall, set S_T_restr_ass
	    and shell_count. */
	
	shell_count = Match_Rsts_Pairs ( s, 
	    obs_restr, S_T_restr_ass,  obs_restr_count, no_restr_count, 
	    pair, S_pair_ass, obs_pair_count, 
	    dist_infl, rev_mapping_table, pair_matched );
	
	if ( shell_count > pair_shell_count[s] ) {
	    nrerror("Error in Process_Distances_Shells: Wrong match");
	}
	
	shell_total_count += shell_count;
	
	/* Calculate the completeness Factor */

	if ( pair_shell_count[s] == 0 ) {
	    compl_fac = NOT_A_NUMBER;
	} else {
	    compl_fac =  shell_count / ( (Distval) pair_shell_count[s] );
	}
	
        if ( terseness < 6 ) {
	    printf("%8d restraints matched to pairs in shell %2d =  %5.1f percent\n\n", 
	            shell_count, s, 100.0 * compl_fac );
        }
	
    }
        
    if ( terseness < 9 ) {
        printf("Found   %8d pairs \n", pair_total_count);
        printf("Found   %8d \"observable\" restraints\n", obs_restr_count - no_restr_count);
        printf("Matched %8d restraints\n", shell_total_count);
        printf("SO      %8d RESTRAINTS ARE NOT OBSERVED IN ENSEMBLE.\n", 
	        obs_restr_count - no_restr_count - shell_total_count );
    }
    
    if ( pair_total_count == 0) {
	total_compl_fac = NOT_A_NUMBER;
    } else {
	total_compl_fac = shell_total_count / (Distval) pair_total_count;
    }
    
    
    ShowSummary ( obs_restr_count, obs_pair_count, S_T_restr_ass, S_pair_ass);\

    ProcessCompletenessOutput(
	dist_count, total_compl_fac, shell_total_count, 
	obs_restr, obs_restr_count, no_restr_count,
	pair, obs_pair_count, no_pair_count,
	S_T_restr_ass, S_pair_ass,
	theSet_ref, theSet, 
	drst_filnam, atominp, obs_atom, obs_atom_count, 
	filnam, dist_infl, rev_mapping_table, rng_filnam );

    		          
    /* Extract pairs that are not restraints */
			    
    /* Mark the pairs that are not matched
     * Afterwards the pairs are rendered useless, don't use them */
     
    no_pair_count += Extract_Pairs ( pair_matched, 
			pair, obs_pair_count, no_pair_count);

    DR_list.list = pair;
    DR_list.dist_count = obs_pair_count;
    DR_list.struc_count = 1;

    /* Write them out again */
    if ( terseness < 6 ) {
        printf("Writing the set of %d not matched pairs.\n", 
	 obs_pair_count-no_pair_count);
    }
    ProcessDROutput( &DR_list, no_pair_count,"noe",false,
	theSet.molecule,theSet.struc_count);


        
    /* Extract restraints from obs_restr that are not matched */
	    
    no_restr_count += Extract_DR ( 
	    S_T_restr_ass, obs_restr, obs_restr_count, no_restr_count );
	    
    DR_list.list = obs_restr;
    DR_list.dist_count = obs_restr_count;
    DR_list.struc_count = 1;

    if ( terseness < 6 ) {
        printf("Writing the set of %d not matched restraints.\n", 
	obs_restr_count-no_restr_count);
    }
    ProcessDROutput( &DR_list, no_restr_count,"noe",false,
	theSet.molecule,theSet.struc_count);
    
    /* After this the list of restraints is rendered useless, don't use it */

    return ( total_compl_fac );
}


static int Extract_Pairs ( Boolean *pair_matched, 
    DRestraint *pair,      int obs_pair_count,  int no_pair_count)
{
    int i, numb_del, no_pair_count_act;
    
    numb_del = 0;
    no_pair_count_act = 0;
        
    for (i=0; i<obs_pair_count; i++) {
    
	if ( !pair[i].keep  ) {            /* Skip "deleted" pairs */
	    no_pair_count_act++;
	    continue;
	}
	
	if ( pair_matched[i] ) {
	    pair[i].keep = false;
	    numb_del++;
	}
    }

    if ( no_pair_count != no_pair_count_act) {
	nrerror("Error in Extract_Pairs: bad no_pair_count_act");
    }

/*    printf ("Will NOT write the %4d matched pairs.\n", numb_del); */

    return ( numb_del );
}



static int Extract_DR ( Shells  *S_T_restr_ass, DRestraint *obs_restr, 
			int obs_restr_count, int no_restr_count )
{
    int i, numb_del, no_restr_count_act;
    int s, t;
    
    numb_del = 0;
    no_restr_count_act = 0;
    
/*
    s = GetOptionIVal( opt_s );
    t = GetOptionIVal( opt_t );
*/
    s = 0;
    t = 0;
    
    if ( ( s < 0 ) || ( t < 0 ) || 
	( s > GetOptionIVal( opt_numb_shells_prs ) ) || 
	( t > GetOptionIVal( opt_numb_shells_rst ) ) ) {
	printf("given shells are s:%d t:%d\n", s, t);
	nrerror("Error in Extract_DR: bad shell info ");
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
	if ( !obs_restr[i].keep  ) {            /* Skip "deleted" restraints */
	    no_restr_count_act++;
	    continue;
	}
	
	if ( ( S_T_restr_ass[i].S != s ) || 
	     ( S_T_restr_ass[i].T != t ) ) {     /* Skip bad categories */
	    obs_restr[i].keep = false;
	    numb_del++;
	}
    }

    if ( no_restr_count != no_restr_count_act) {
	nrerror("Error in Partition_Restraints: bad no_restr_count_act");
    }

    if ( terseness < 6 ) {
        printf ("Will write %4d restraints matched\n", 
	obs_restr_count-(no_restr_count+numb_del));
    }
    return ( numb_del );
}



static void ProcessCompletenessOutput(
	int dist_count, Distval total_compl_fac, int shell_total_count, 
	DRestraint  *obs_restr, int obs_restr_count, int no_restr_count,
	DRestraint  *pair,      int obs_pair_count,  int no_pair_count,
	Shells *S_T_restr_ass, int *S_pair_ass,
	Strucset theSet_ref, Strucset theSet, 
	Filnam drst_filnam, char *atominp, 
	Boolean *obs_atom, int obs_atom_count, 
	Filnam filnam, Distinfl *dist_infl, int *rev_mapping_table, Filnam rng_filnam )
{

    FILE     *fp;
    Qtyp     qto;     /* NEW !! */
    Distval avCc, avCa, avCr, avCaT, avCrT, avCar;
    Distval sdCc, sdCa, sdCr, sdCaT, sdCrT, sdCar;
    
    avCc=0, avCa=0, avCr=0, avCaT=0, avCrT=0, avCar=0;
    sdCc=0, sdCa=0, sdCr=0, sdCaT=0, sdCrT=0, sdCar=0;
    
    
/* open the completeness output file using the ident of a previous file */

    if ( terseness < 6 ) {
        puts( "Creating the restraint completeness output files" );
    }
    
    qto = noe_compl_data;
    fp = OpenAquaFile( qto, "w", filnam, true );
    if ( !fp )
        exit( 1 );
    WriteAquaIDLines( fp, qto, filnam );
    fprintf( fp, "\n$ Restraints were taken from: %s\n", drst_filnam );
    WriteTxtBuffer( fp );
    

    Show_Explain_AQUA ( fp );
    Show_General_AQUA ( fp, dist_count, theSet );

    Show_Paramet_AQUA ( fp, atominp, rng_filnam );

    Show_OverSum_AQUA ( fp, total_compl_fac, 
	    obs_pair_count,  no_pair_count,
	    obs_restr_count, no_restr_count, shell_total_count, 
	    obs_atom_count, theSet );
    
    Show_Summary_AQUA ( fp, 
		    obs_restr_count, obs_pair_count, 
		    S_T_restr_ass, S_pair_ass );
    
    Show_PerClass_AQUA ( fp,
	    obs_restr, obs_restr_count, no_restr_count,
	    pair,      obs_pair_count,  no_pair_count,
	    S_T_restr_ass, 
	    &avCc, &sdCc );
	 
    Show_PerAtom_Type_AQUA ( fp,
	    obs_restr, obs_restr_count, no_restr_count,
	    pair,      obs_pair_count,  no_pair_count,
	    S_T_restr_ass, 
	    &avCaT, &sdCaT );
	    
    Show_PerResidue_Type_AQUA ( fp,
	    obs_restr, obs_restr_count, no_restr_count,
	    pair,      obs_pair_count,  no_pair_count,
	    theSet, theSet_ref, S_T_restr_ass, 
	    &avCrT, &sdCrT );

    Show_PerResidue_AQUA ( fp,
	    obs_restr, obs_restr_count, no_restr_count,
	    pair,      obs_pair_count,  no_pair_count,
	    theSet, theSet_ref, S_T_restr_ass, 
	    &avCr, &sdCr );


    if ( OptionTrue( opt_detail ) ) {

	Show_PerResAtom_AQUA ( fp,
		obs_restr, obs_restr_count, no_restr_count,
		pair,      obs_pair_count,  no_pair_count,
		theSet, theSet_ref, S_T_restr_ass,
		dist_infl, rev_mapping_table, dist_count,
		&avCar, &sdCar );
		
	Show_PerAtom_AQUA ( fp,
		obs_restr, obs_restr_count, no_restr_count,
		pair,      obs_pair_count,  no_pair_count,
		theSet, S_T_restr_ass, obs_atom, obs_atom_count, 
		dist_infl, rev_mapping_table, dist_count,
		&avCa, &sdCa );
    }
           

    Show_AVs_SDs_AQUA ( fp,
	avCc, avCaT, avCrT, avCa, avCr, avCar, 
	sdCc, sdCaT, sdCrT, sdCa, sdCr, sdCar );

   fclose( fp );
}


static void Show_Explain_AQUA( FILE *fp ) 
{
    fputc( '\n', fp );
    fprintf( fp, "$ EXPLAIN\n" );

    fprintf( fp, "$ Pairs are theoretically observable, as based on the structure(s)\n");
    fprintf( fp, "$ whereas restraints have experimentally been observed.\n" );



    fprintf( fp, "$ COMBINATIONS of the following explained:\n" );
    fprintf( fp, "$ C      : completeness           cum    : cumulative\n" );
    fprintf( fp, "$ M      : matched to pairs       av     : average\n" );
    fprintf( fp, "$ N      : number of              sd     : standard deviation\n" );
    fprintf( fp, "$ O      : observable             n      : not matched\n" );
    fprintf( fp, "$ P      : pairs                  Low    : lower distance\n" );
    fprintf( fp, "$ R      : restraints             Upp    : upper distance\n" );
    fprintf( fp, "$ S      : shells\n" );
    fprintf( fp, "$ \n" );
    
    fprintf( fp, "$ END\n" );
}

static void Show_General_AQUA( FILE *fp, int dist_count, Strucset theSet ) 
{
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  GENERAL\n" );
    fprintf( fp, "$ LIST #RECORDS %i\n", 3 );
    fprintf( fp, "NRestraints %i\n", dist_count );
    fprintf( fp, "NStructures %i\n", theSet.struc_count );
    fprintf( fp, "NResidues   %i\n", theSet.molecule->res_count );
    fprintf( fp, "$ END\n" );
}



static void Show_Paramet_AQUA( FILE *fp, Filnam filnam_ob, Filnam rng_filnam ) 
{
    fputc( '\n', fp );
    fprintf( fp, "$ DATA PARAMETERS\n" );
    fprintf( fp, "$ LIST #RECORDS %i\n", 13 );
    
    fprintf( fp, "Low_P                         %.2f\n",  GetOptionRVal( opt_min_dist_prs ) );
    fprintf( fp, "Upp_P                         %.2f\n",  GetOptionRVal( opt_max_dist_prs ) );
    fprintf( fp, "NS_P                          %d\n",    GetOptionIVal( opt_numb_shells_prs ) );

    fprintf( fp, "Low_R                         %.2f\n",  GetOptionRVal( opt_min_dist_rst ) );
    fprintf( fp, "Upp_R                         %.2f\n",  GetOptionRVal( opt_max_dist_rst ) );
    fprintf( fp, "NS_R                          %d\n",    GetOptionIVal( opt_numb_shells_rst ) );

    fprintf( fp, "observable_atoms_file         %s\n",    filnam_ob );
    fprintf( fp, "range_file                    %s\n",    rng_filnam );
    
    fprintf( fp, "use_intra-residuals           %c\n",    Boolean_To_Char( OptionTrue( opt_use_intra ) ) );
    fprintf( fp, "threshold_redundancy          %.0f\n",  GetOptionRVal( opt_thresh ) );
    fprintf( fp, "averaging power               %.2f\n",  GetOptionRVal( opt_avg ) );
    fprintf( fp, "savpseu                       %c\n",    Boolean_To_Char( OptionTrue( opt_savpseu ) ) );
    fprintf( fp, "double_count                  %c\n",    Boolean_To_Char( OptionTrue( opt_double_count ) ) );
    fprintf( fp, "analyse_ssa                   %c\n",    Boolean_To_Char( OptionTrue( opt_anal_ssa ) ) );
    
    fprintf( fp, "$ END\n" );
}


static void Show_OverSum_AQUA( FILE *fp, Distval total_compl_fac, 
	int obs_pair_count,  int no_pair_count,
	int obs_restr_count, int no_restr_count, int shell_total_count,
	int obs_atom_count,  Strucset theSet) 
{
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  SUMMARY_OVERALL\n" );
    fprintf( fp, "$ REMARK O_stru : number of observable atoms in test structure\n" );
    fprintf( fp, "$ REMARK A_stru : number of all atoms in test structure (including pseudos)\n" );
    fprintf( fp, "$ TABLE #RECORDS 1 #ITEMS 6\n");
    fprintf( fp, "     C    N_P   N_RM  N_RnM    N_R O_stru A_stru\n");
    fprintf( fp, "%6.1f %6i %6i %6i %6i %6i %6i\n",
		 100.0*total_compl_fac, 
		 obs_pair_count  - no_pair_count, 
		 shell_total_count,
		 obs_restr_count - no_restr_count - shell_total_count,
		 obs_restr_count - no_restr_count,
		 obs_atom_count, theSet.molecule->atom_count );
    fprintf( fp, "$ END\n" );
}



static void Show_AVs_SDs_AQUA( FILE *fp, 
    Distval avCc, Distval avCaT, Distval avCrT, Distval avCa, Distval avCr, Distval avCra, 
    Distval sdCc, Distval sdCaT, Distval sdCrT, Distval sdCa, Distval sdCr, Distval sdCra )
{
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  COMPLETENESS_AVERAGES_AND_SDS\n" );
    fprintf( fp, "$ a      : on atom basis          r      : on residue basis\n" );
    fprintf( fp, "$ c      : on class basis         t      : on type basis\n" );    
    fprintf( fp, "$ TABLE #RECORDS 2 #ITEMS 7\n");
    fprintf( fp, "Quantity       c     at     rt      r     ra      a\n");
    fprintf( fp, "av        %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f\n", 
		    avCc, avCaT, avCrT, avCr, avCra, avCa);
    fprintf( fp, "sd        %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f\n", 
		    sdCc, sdCaT, sdCrT, sdCr, sdCra, sdCa );
    fprintf( fp, "$ END\n" );
}



static void Show_PerClass_AQUA ( FILE *fp,
	DRestraint  *obs_restr, int obs_restr_count, int no_restr_count,
	DRestraint  *pair,      int obs_pair_count,  int no_pair_count,
	Shells *S_T_restr_ass,  
	Distval *av, Distval *sd )
{
			       /* See Sort_Restraints routine */
#define NUM_LIST      5

    int i, skip_count;
    int *rst_count, *no_rst_match_count, *prs_count, *dummy;
    char classes[ NUM_LIST  ][ 11 ] =
	{"intra", "sequential", "medium", "long", "inter"};
    Shells *fake = NULL;
    Distval *completeness, Nsd;
    
    skip_count = 0;
    completeness = CREATE_NEW( Distval, NUM_LIST);
    
    Partition_Restraints_Classes ( obs_restr, obs_restr_count, no_restr_count, 
	    &rst_count, true,  &no_rst_match_count, S_T_restr_ass );
    Partition_Restraints_Classes ( pair,      obs_pair_count,  no_pair_count,  
	    &prs_count, false, &dummy,              fake );

    for (i=0; i<NUM_LIST; i++) {
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] = 100.0 * rst_count[i] / prs_count[i];
	}	
    }
    
    *av = Average_Float_Array( completeness, NUM_LIST, skip_count );    
    *sd = SD_Float_Array( completeness, NUM_LIST, skip_count );
			    
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  COMPLETENESS_PER_CLASS\n" );
    fprintf( fp, "$ TABLE #RECORDS 5 #ITEMS 7\n");
    fprintf( fp, "Type          N_P   N_RM      C  N_RnM  N_sd Remark\n");
    
    for (i=0; i<NUM_LIST; i++) {
    
        fprintf( fp, "%-10s ", classes[i] );
	fprintf( fp, "%6i ", prs_count[i] );
	fprintf( fp, "%6i ", rst_count[i] );
	fprintf( fp, "%6.1f ", completeness[i] );
	fprintf( fp, "%6i ", no_rst_match_count[i] );

	if ( (Same_Floats (completeness[i], (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) 0.0 )) ) {
	    Nsd = NOT_A_NUMBER;
	} else {
	    Nsd = ( completeness[i] - *av ) / *sd;
	}
	fprintf( fp, "%5.1f ",  Nsd);

	if ( ( i==0 ) && !OptionTrue( opt_use_intra ) ) {
	    fprintf( fp, "Not_analysed\n" );
	} else {
	    if ( ( i==NUM_LIST-1 ) && ( prs_count[i] == 0 ) ) {
		fprintf( fp, "No_chains\n" );
	    } else {
		if ( ( Nsd < -1.0 ) && 
		     ( ! Same_Floats (Nsd, (Distval) NOT_A_NUMBER) ) ) {
		    fprintf( fp, ">sigma\n" );
		} else {
		    fprintf( fp, "_\n" );
		}
	    }
	}
    }
    fprintf( fp, "$ END\n" );
        
    free ( rst_count );
    free ( no_rst_match_count );
    free ( prs_count );
    free ( completeness );
#undef NUM_LIST
}	



static void Show_PerAtom_Type_AQUA ( FILE *fp,
	DRestraint  *obs_restr, int obs_restr_count, int no_restr_count,
	DRestraint  *pair,      int obs_pair_count,  int no_pair_count,
	Shells *S_T_restr_ass, Distval *av, Distval *sd )
{

#define NUM_LIST 8

    int i, skip_count;
    int *rst_count, *no_rst_match_count, *prs_count, *dummy;
    Shells *fake = NULL;
    Distval *completeness, Nsd;

    char classes[ NUM_LIST ][ 8 ] = 
	{ "amide", "alpha", "beta", "gamma",
	  "delta", "epsilon", "dzeta  ", "other"};
	  
    
    skip_count = 0;
    completeness = CREATE_NEW( Distval, NUM_LIST);
    
    Partition_Restraints_AtomT ( obs_restr, obs_restr_count, no_restr_count, 
	    &rst_count, true,  &no_rst_match_count, S_T_restr_ass );
    Partition_Restraints_AtomT ( pair,      obs_pair_count,  no_pair_count,  
	    &prs_count, false, &dummy,              fake );
			    
    for (i=0; i<NUM_LIST; i++) {
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] =  (Distval) 100.0 * rst_count[i] / 
				    prs_count[i];
	}	
    }
    
    *av = Average_Float_Array( completeness, NUM_LIST, skip_count );    
    *sd = SD_Float_Array( completeness, NUM_LIST, skip_count );

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  COMPLETENESS_PER_ATOM_TYPE\n" );
    fprintf( fp, "$ TABLE #RECORDS 8 #ITEMS 7\n");
    fprintf( fp, "Type          N_P   N_RM      C  N_RnM  N_sd Remark\n");
    
    for (i=0; i<NUM_LIST; i++) {
    
        fprintf( fp, "%-10s ", classes[i] );
	fprintf( fp, "%6i ", prs_count[i] );
	fprintf( fp, "%6i ", rst_count[i] );
	fprintf( fp, "%6.1f ", completeness[i] );
	fprintf( fp, "%6i ", no_rst_match_count[i] );
		
	if ( (Same_Floats (completeness[i], (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) 0.0 )) ) {
	    Nsd = NOT_A_NUMBER;
	} else {
	    Nsd = ( completeness[i] - *av ) / *sd;
	}
	fprintf( fp, "%5.1f ",  Nsd);
	if (( Nsd < -1.0 ) && ( ! Same_Floats (Nsd, (Distval) NOT_A_NUMBER) ) ) {
	    fprintf( fp, ">sigma\n" );
	} else {
	    fprintf( fp, "_\n" );
	}
    }
    
    fprintf( fp, "$ END\n" );
        
    free ( rst_count );
    free ( no_rst_match_count );
    free ( prs_count );
    free ( completeness );

#undef NUM_LIST
}	



static void Show_PerResidue_Type_AQUA ( FILE *fp,
	DRestraint  *obs_restr, int obs_restr_count, int no_restr_count,
	DRestraint  *pair,      int obs_pair_count,  int no_pair_count,
	Strucset theSet, Strucset theSet_ref, Shells *S_T_restr_ass,
	Distval *av, Distval *sd )
{

    int i, j, skip_count, num_list;
    int *rst_count, *no_rst_match_count, *prs_count, *dummy;
    int *res_count, *ato_count;
    Shells *fake = NULL;
    Distval *completeness, Nsd;
    ID  *resID;
    
    skip_count = 0;
    num_list = theSet_ref.molecule->res_count;
    
    resID = CREATE_NEW( ID, num_list);
    for (i=0; i<num_list; i++) {
	memcpy( resID[i], theSet_ref.molecule->residues[i].res_nam, size_ID );
    }
    
    completeness = CREATE_NEW( Distval, num_list);
    res_count = CREATE_NEW( int, num_list);
    ato_count = CREATE_NEW( int, num_list);
    
    Partition_Restraints_ResidueT ( obs_restr, obs_restr_count, no_restr_count, 
	    &rst_count, true,  &no_rst_match_count, S_T_restr_ass, 
	    resID, num_list );
    Partition_Restraints_ResidueT ( pair,      obs_pair_count,  no_pair_count,  
	    &prs_count, false, &dummy,              fake, 
	    resID, num_list );

    for (i=0; i<num_list; i++) {
    	
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] = 100.0 * rst_count[i] / prs_count[i];
	}

	res_count[i] = 0;
	for (j=0; j< theSet.molecule->res_count; j++) {
	    if ( EQID ( theSet.molecule->residues[j].res_nam, resID[ i ] ))  {
		res_count[i]++;
	    }
	}

	ato_count[i] = 0;	
	for (j=0; j< theSet_ref.molecule->atom_count; j++) {
	    if ( EQID ( theSet_ref.molecule->atoms[j].res_nam, resID[i] ) ) {
		ato_count[i]++;
	    }
	}
    }
    
    *av = Average_Float_Array( completeness, num_list, skip_count );    
    *sd = SD_Float_Array( completeness, num_list, skip_count );
			    
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  COMPLETENESS_PER_RESIDUE_TYPE\n" );
    fprintf( fp, "$ TABLE #RECORDS %d #ITEMS 9\n", num_list);
    fprintf( fp, "Resnm    N_r  N_O_a    N_P   N_RM      C  N_RnM  N_sd Remark\n");
    
    for (i=0; i<num_list; i++) {
    
        FPRINTID(fp, resID[i]);
	fprintf( fp, "   ");
	fprintf( fp, "%6i ", res_count[i] );
	fprintf( fp, "%6i ", ato_count[i] );
	fprintf( fp, "%6i ", prs_count[i] );
	fprintf( fp, "%6i ", rst_count[i] );
	fprintf( fp, "%6.1f ", completeness[i] );
	fprintf( fp, "%6i ", no_rst_match_count[i] );
		
	if ( (Same_Floats (completeness[i], (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) 0.0 )) ) {
	    Nsd = NOT_A_NUMBER;
	} else {
	    Nsd = ( completeness[i] - *av ) / *sd;
	}
	fprintf( fp, "%5.1f ",  Nsd);
	if (( Nsd < -1.0 ) && ( ! Same_Floats (Nsd, (Distval) NOT_A_NUMBER) ) ) {
	    fprintf( fp, ">sigma\n" );
	} else {
	    fprintf( fp, "_\n" );
	}
    }
    
    fprintf( fp, "$ END\n" );
        
    free ( rst_count );
    free ( no_rst_match_count );
    free ( prs_count );
    free ( completeness );
    free ( res_count );
    free ( ato_count );
    free ( resID );

}	


static void Show_PerResidue_AQUA ( FILE *fp,
	DRestraint  *obs_restr, int obs_restr_count, int no_restr_count,
	DRestraint  *pair,      int obs_pair_count,  int no_pair_count,
	Strucset theSet, Strucset theSet_ref, Shells *S_T_restr_ass, 
	Distval *av, Distval *sd )
{

    int i, j, skip_count, num_list;
    int *rst_count, *no_rst_match_count, *prs_count, *dummy, *ato_count;
    Shells *fake = NULL;
    Distval *completeness, Nsd;
    
    skip_count = 0;
    num_list = theSet.molecule->res_count;
    completeness = CREATE_NEW( Distval, num_list);
    ato_count = CREATE_NEW( int, num_list);
    
    Partition_Restraints_Residue ( obs_restr, obs_restr_count, no_restr_count, 
	    &rst_count, true,  &no_rst_match_count, S_T_restr_ass, theSet );
    Partition_Restraints_Residue ( pair,      obs_pair_count,  no_pair_count,  
	    &prs_count, false, &dummy,              fake, theSet );
			    
    for (i=0; i<num_list; i++) {
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] = 100.0 * rst_count[i] / prs_count[i];
	}	

	ato_count[i] = 0;	
	for (j=0; j< theSet_ref.molecule->atom_count; j++) {
	    if ( EQID ( theSet_ref.molecule->atoms[j].res_nam, 
	                 theSet.molecule->residues[i].res_nam ) ) {
		ato_count[i]++;
	    }
	}

    }
    
    *av = Average_Float_Array( completeness, num_list, skip_count );    
    *sd = SD_Float_Array( completeness, num_list, skip_count );

    fputc( '\n', fp );
    fprintf( fp, "$ DATA  COMPLETENESS_PER_RESIDUE\n" );
    fprintf( fp, "$ TABLE #RECORDS %d #ITEMS 10\n", num_list);
    fprintf( fp, "Chn Resnm Res  N_O_a    N_P   N_RM      C  N_RnM  N_sd Remark\n");	

    for (i=0; i<num_list; i++) {
        WriteFResidue ( fp, &theSet.molecule->residues[ i ] ); /* 3 items */
	fprintf( fp, "%6i ", ato_count[i] );
	fprintf( fp, "%6i ", prs_count[i] );
	fprintf( fp, "%6i ", rst_count[i] );
	
	fprintf( fp, "%6.1f ", completeness[i] );
	
	fprintf( fp, "%6i ", no_rst_match_count[i] );

	if ( (Same_Floats (completeness[i], (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) 0.0 )) ) {
	    Nsd = NOT_A_NUMBER;
	} else {
	    Nsd = ( completeness[i] - *av ) / *sd;
	}
	fprintf( fp, "%5.1f ",  Nsd);
	if (( Nsd < -1.0 ) && ( ! Same_Floats (Nsd, (Distval) NOT_A_NUMBER) ) ) {
	    fprintf( fp, ">sigma\n" );
	} else {
	    fprintf( fp, "_\n" );
	}
    }
    
    fprintf( fp, "$ END\n" );
        
    free ( rst_count );
    free ( no_rst_match_count );
    free ( prs_count );
    free ( completeness );
    free ( ato_count );
}	



static void Show_PerAtom_AQUA ( FILE *fp,
	DRestraint  *obs_restr, int obs_restr_count, int no_restr_count,
	DRestraint  *pair,      int obs_pair_count,  int no_pair_count,
	Strucset theSet, Shells *S_T_restr_ass,
	Boolean *obs_atom, int obs_atom_count, 
	Distinfl *dist_infl, int *rev_mapping_table, int dist_count,
	Distval *av, Distval *sd )
{

    int	    i, k, skip_count, *list_idx;
    int	    *rst_count, *no_rst_match_count, *prs_count, *dummy;
    Shells  *fake = NULL;
    Distval   *completeness, Nsd;
    
    skip_count = 0;
    completeness = CREATE_NEW( Distval, obs_atom_count);
    list_idx     = CREATE_NEW(   int, obs_atom_count );
    
    k=0;
    
    for (i=0; i<theSet.molecule->atom_count; i++) {
	if ( obs_atom[i] ) {
	    list_idx[ k++ ] = i;
	}
    }
    
    if ( k != obs_atom_count ) {
	nrerror("Error in Show_PerAtom_AQUA: bad count ");
    }

    Partition_Restraints_Atom ( obs_restr, obs_restr_count, no_restr_count, 
	    &rst_count, true,  &no_rst_match_count, S_T_restr_ass, list_idx,
	    obs_atom_count, dist_infl, rev_mapping_table, dist_count);
	    
    Partition_Restraints_Atom ( pair,      obs_pair_count,  no_pair_count,
	    &prs_count, false, &dummy,              fake, list_idx,
	    obs_atom_count, dist_infl, rev_mapping_table, dist_count );

    for (i=0; i<obs_atom_count; i++) {
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] = 100.0 * rst_count[i] / prs_count[i];
	}	
    }
    
    *av = Average_Float_Array( completeness, obs_atom_count, skip_count );    
    *sd = SD_Float_Array( completeness, obs_atom_count, skip_count );
			    
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  COMPLETENESS_PER_ATOM\n" );
    fprintf( fp, "$ TABLE #RECORDS %d #ITEMS 7\n", obs_atom_count);
    fprintf( fp, "Chn Resnm Res Atom    N_P   N_RM      C  N_RnM  N_sd Remark\n");
    
    for (i=0; i<obs_atom_count; i++) {
        
	WriteFAtom( fp, &theSet.molecule->atoms[ list_idx[i] ], use_int );
	
	fprintf(   fp, "%6i ", prs_count[i] );
	fprintf(   fp, "%6i ", rst_count[i] );
	
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] = 100.0 * rst_count[i] / prs_count[i];
	}
	fprintf( fp, "%6.1f ", completeness[i] );
	fprintf( fp, "%6i ", no_rst_match_count[i] );

	if ( (Same_Floats (completeness[i], (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) 0.0 )) ) {
	    Nsd = NOT_A_NUMBER;
	} else {
	    Nsd = ( completeness[i] - *av ) / *sd;
	}
	fprintf( fp, "%5.1f ",  Nsd);
	if (( Nsd < -2.0 ) && ( ! Same_Floats (Nsd, (Distval) NOT_A_NUMBER) ) ) {
	    fprintf( fp, ">2sigma\n" );
	} else {
	    fprintf( fp, "_\n" );
	}
    }
    
    fprintf( fp, "$ END\n" );
        
    free ( rst_count );
    free ( no_rst_match_count );
    free ( prs_count );
    free ( completeness );
    free ( list_idx );
}	



static void Show_PerResAtom_AQUA ( FILE *fp,
	DRestraint  *obs_restr, int obs_restr_count, int no_restr_count,
	DRestraint  *pair,      int obs_pair_count,  int no_pair_count,
	Strucset theSet, Strucset theSet_ref, Shells *S_T_restr_ass,
	Distinfl *dist_infl, int *rev_mapping_table, int dist_count,
	Distval *av, Distval *sd )
{
/* Show the data per atom for each residue once ! */
    int	    i, j, skip_count, *ato_count;
    int	    *rst_count, *no_rst_match_count, *prs_count, *dummy;
    int     obs_atom_count_ref;
    Shells  *fake = NULL;
    Distval   *completeness, Nsd;
    
    
    obs_atom_count_ref = theSet_ref.molecule->atom_count;
    
    skip_count = 0;
    completeness = CREATE_NEW( Distval, obs_atom_count_ref);
    ato_count = CREATE_NEW( int, obs_atom_count_ref);

    Partition_Restraints_ResAtom ( obs_restr, obs_restr_count, no_restr_count, 
	    &rst_count, true,  &no_rst_match_count, S_T_restr_ass,
	    theSet, theSet_ref, dist_infl, rev_mapping_table, dist_count);
	    
    Partition_Restraints_ResAtom ( pair,      obs_pair_count,  no_pair_count,
	    &prs_count, false, &dummy,              fake,
	    theSet, theSet_ref, dist_infl, rev_mapping_table, dist_count );

    for (i=0; i<obs_atom_count_ref; i++) {
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] = 100.0 * rst_count[i] / prs_count[i];
	}	
    }
    
    *av = Average_Float_Array( completeness, obs_atom_count_ref, skip_count );    
    *sd = SD_Float_Array( completeness, obs_atom_count_ref, skip_count );
			    
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  COMPLETENESS_PER_RESATOM\n" );
    fprintf( fp, "$ TABLE #RECORDS %d #ITEMS 11\n", obs_atom_count_ref);
    fprintf( fp, "Chn Resnm Res Atom  N_O_a    N_P   N_RM      C  N_RnM  N_sd Remark\n");
    
    for (i=0; i<obs_atom_count_ref; i++) {
        
	WriteFAtom( fp, &theSet_ref.molecule->atoms[ i ], use_int );
	
	ato_count[i] = 0;	
	for (j=0; j< theSet.molecule->res_count; j++) {
	    if ( EQID ( theSet_ref.molecule->atoms[ i ].res_nam,
	                theSet.molecule->residues[  j ].res_nam ) ) {
		ato_count[i]++;
	    }
	}
	
	fprintf(   fp, "%6i ", ato_count[i] );
	
	fprintf(   fp, "%6i ", prs_count[i] );
	fprintf(   fp, "%6i ", rst_count[i] );
	
	if ( prs_count[i] == 0 ) {
	    completeness[i] = NOT_A_NUMBER;
	    skip_count++;
		     /* Allows a calculation of average and s.d.s */
	} else {
	    completeness[i] = 100.0 * rst_count[i] / prs_count[i];
	}
	fprintf( fp, "%6.1f ", completeness[i] );
	fprintf( fp, "%6i ", no_rst_match_count[i] );

	if ( (Same_Floats (completeness[i], (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) NOT_A_NUMBER)) || 
	     (Same_Floats (*sd, (Distval) 0.0 )) ) {
	    Nsd = NOT_A_NUMBER;
	} else {
	    Nsd = ( completeness[i] - *av ) / *sd;
	}
	fprintf( fp, "%5.1f ",  Nsd);
	if (( Nsd < -1.0 ) && ( ! Same_Floats (Nsd, (Distval) NOT_A_NUMBER) ) ) {
	    fprintf( fp, ">sigma\n" );
	} else {
	    fprintf( fp, "_\n" );
	}
    }
    
    fprintf( fp, "$ END\n" );
        
    free ( rst_count );
    free ( no_rst_match_count );
    free ( prs_count );
    free ( completeness );
    free ( ato_count );
}	



static Boolean Same_Floats (Distval a, Distval b) 
{
    Distval temp;
    
    if ( a > b ) {
	temp = a - b;
    } else {
	temp = b - a;
    }
    
    if ( temp < SMALL_POS_DOUBLE_NUMBER ) {
	return ( true );
    } else {
	return ( false );	
    }
}


static Distval Average_Float_Array ( Distval *f_array, 
	int count, int skip_count )
{
    int i, skip_count_actual;
    Distval result, sum;
    
    result            = 0.0;
    sum               = 0.0;
    skip_count_actual = 0;
   
    if ( count == skip_count ) {
	if ( OptionTrue( opt_detail ) ) {
	    printf ("Warning in Average_Float_Array: Nothing to average\n");
	}
	return ( NOT_A_NUMBER );
    }
    
    for (i=0; i<count; i++) {
	if ( Same_Floats (f_array[i], (Distval) NOT_A_NUMBER) ) {
	    skip_count_actual++;
	    continue;
	}
	sum += f_array[i];
    }

    if ( skip_count_actual != skip_count ) {
	nrerror ("Wrong code");    
    }

    result = sum / ( count - skip_count );

    return ( result );
}



static Distval SD_Float_Array ( Distval *f_array, 
	int count, int skip_count )
{
    int i, skip_count_actual;
    Distval sum, result, temp;
    Distval average = NOT_A_NUMBER; /* To prevent compiler warning */
    
    result            = 0.0;
    sum               = 0.0;
    skip_count_actual = 0;
    
    if ( count == skip_count ) {
	if ( OptionTrue( opt_detail ) ) {
	    printf ("Warning in SD_Float_Array: Nothing to average\n");
	}
	return ( NOT_A_NUMBER );
    }

    for (i=0; i<count; i++) {
	if ( Same_Floats (f_array[i], (Distval) NOT_A_NUMBER) ) {
	    skip_count_actual++;
	    continue;
	}
	sum += f_array[i];
    }

    if ( skip_count_actual != skip_count ) {
	nrerror ("Wrong code");    
    }

    if ( ( count - skip_count ) == 0 ) {
	nrerror ("Error in Average_Float_Array: Nothing to average");
    } else {
	average = sum / ( count - skip_count );
    }
    
    sum = 0;
    
    for (i=0; i<count; i++) {
	if ( f_array[i] < 0.0 ) {
	    continue;
	}
        temp =  average - f_array[i];
	sum += temp * temp;
    }

    result = sqrt ( sum / ( count - skip_count ) ) ;
    
    return ( result );
}




static void Partition_Restraints_Classes ( 
    DRestraint *obs_restr, int obs_restr_count, int no_restr_count,
    int **rst_count, Boolean do_no_rst_count, int **no_rst_match_count, 
    Shells *S_T_restr_ass )
{

/* See sorting code in Sort_Restraints routine */
#define MEDIUM_BORDER 4
#define NUM_LIST      5

    typedef enum { intra, seq, medium, llong, inter } Contact_Range;

    char    *fct_name = "Partition_Restraints_Classes";
    int		    i, no_restr_count_act, diff_resnum;
    Contact_Range   contact_range;
    Boolean	    matched;
    
    *rst_count          = CREATE_NEW( int, NUM_LIST );
    *no_rst_match_count = CREATE_NEW( int, NUM_LIST );
    
    no_restr_count_act = 0;
    	
    for (i=0; i<NUM_LIST; i++) {
	(*rst_count)[i] = 0;
	(*no_rst_match_count)[i] = 0;
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
        matched = false;

        if ( !obs_restr[i].keep ) {
	    no_restr_count_act++;
	    continue;
	}

	if ( do_no_rst_count ) {
	    if ( S_T_restr_ass[i].S == 0 ) {
	    
		if ( S_T_restr_ass[i].T != 0 ) {	    
		    printf("restraint %d has flags S and T; %d %d\n", i, 
			    S_T_restr_ass[i].S, S_T_restr_ass[i].T);
		    nrerror("Error in Partition_Restraints: bad info ");
		}
		matched = false;
		
	    } else {
	    	matched = true;
	    }
	}
	
	diff_resnum = obs_restr[i].atom_A.res_num_orig - 
	              obs_restr[i].atom_B.res_num_orig ;
	
	if (diff_resnum < 0)
	    diff_resnum *= -1; 
	   	
	contact_range = intra; /* Default */
	    
	if ( diff_resnum == 1 )
	    contact_range = seq;
	    
	if ( ( diff_resnum > 1 ) && ( diff_resnum <= MEDIUM_BORDER ) )
	    contact_range = medium;
	    
	if ( diff_resnum > MEDIUM_BORDER )
	    contact_range = llong;
	    
	if ( NOTEQID( obs_restr[i].atom_A.chain, obs_restr[i].atom_B.chain ) )
	    contact_range = inter;

        if ( ( contact_range == intra ) && ( diff_resnum != 0 ) ) {
            FATAL("Code error");
        }

	switch ( contact_range ) {
	
	    case intra :
		if ( do_no_rst_count ) {
		    if ( !matched ) {
			(*no_rst_match_count)[0]++;
		    } else {
			++(*rst_count)[0];
		    }
		} else {
		    ++(*rst_count)[0];
		}
		break;
	    case seq :
		if ( do_no_rst_count ) {
		    if ( !matched ) {
			(*no_rst_match_count)[1]++;
		    } else {
			++(*rst_count)[1];
		    }
		} else {
		    ++(*rst_count)[1];
		}
		break;
	    case medium :
		if ( do_no_rst_count ) {
		    if ( !matched ) {
			(*no_rst_match_count)[2]++;
		    } else {
			++(*rst_count)[2];
		    }
		} else {
		    ++(*rst_count)[2];
		}
		break;
	    case llong :
		if ( do_no_rst_count ) {
		    if ( !matched ) {
			(*no_rst_match_count)[3]++;
		    } else {
			++(*rst_count)[3];
		    }
		} else {
		    ++(*rst_count)[3];
		}
		break;
	    case inter :
		if ( do_no_rst_count ) {
		    if ( !matched ) {
			(*no_rst_match_count)[4]++;
		    } else {
			++(*rst_count)[4];
		    }
		} else {
		    ++(*rst_count)[4];
		}
		break;
	    default :
		nrerror("ERROR in Partition_Restraints: default case");
	}

    }
    
    if ( no_restr_count != no_restr_count_act) {
	nrerror("Error in Partition_Restraints: bad no_restr_count_act");
    }
#undef NUM_LIST
#undef MEDIUM_BORDER
}			   


static void Partition_Restraints_AtomT ( 
    DRestraint *obs_restr, int obs_restr_count, int no_restr_count,
    int **rst_count, Boolean do_no_rst_count, int **no_rst_match_count, 
    Shells *S_T_restr_ass )
{

#define NUM_LIST 8

    char AtomChar[ NUM_LIST ] = { '_', 'A', 'B', 'G', 'D', 'E', 'Z', ' '};
    
    int		i, no_restr_count_act, atom_idx;
    int		list_number, atom_a_list_number;
    Boolean	double_count, matched;
    ID		id;
    
    *rst_count          = CREATE_NEW( int, NUM_LIST );
    *no_rst_match_count = CREATE_NEW( int, NUM_LIST );
    
    double_count = OptionTrue( opt_double_count );
    
    no_restr_count_act = 0;
    	
    for (i=0; i<NUM_LIST; i++) {
	(*rst_count)[i] = 0;
	(*no_rst_match_count)[i] = 0;
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
        matched = false;

        	if ( !obs_restr[i].keep ) {
	    no_restr_count_act++;
	    continue;
	}

	if ( do_no_rst_count ) {
	    if ( S_T_restr_ass[i].S == 0 ) {
	    
		if ( S_T_restr_ass[i].T != 0 ) {	    
		    printf("restraint %d has flags S and T; %d %d\n", i, 
			    S_T_restr_ass[i].S, S_T_restr_ass[i].T);
		    nrerror("Error in Partition_Restraints_AtomT: bad info ");
		}
		matched = false;
		
	    } else {
	    	matched = true;
	    }
	}
	
	for ( atom_idx = 0; atom_idx < 2; atom_idx++ ) {
	
	    if ( atom_idx == 0 ) {
	        CopyChars( id, obs_restr[i].atom_A.atom_int, size_ID );
	    } else {
	        CopyChars( id, obs_restr[i].atom_B.atom_int, size_ID );
	    }
	    
        /* Special casing the amide proton */
        if ( ( LenChars( id, size_ID ) == 1 ) && ( id[0] == 'H' ) ) {
            list_number = 0;
/*
            printf("%c  | %c\n", id[0], AtomChar[list_number]);
*/
        } else {
	        for ( list_number = 1; list_number < (NUM_LIST-1); list_number++ ) {
		       if ( memchr ( id, AtomChar[ list_number ], size_ID ) ) {                   
/*                   
                   printf("%c%c | %c\n", id[0], id[1], AtomChar[list_number]);
*/
		           break;
		       }
	        }
        }
	    
	    /* Prevent a double count in same list member ? */
	                    
            atom_a_list_number = -1; /* default but just inserted for compiler warnings */
            
	    if ( double_count ) {
		if ( atom_idx == 0 ) {
		    atom_a_list_number = list_number;
		} else { 
		    if ( atom_a_list_number == list_number) {
			break;
		    }
		}	    
	    }
	    
	    if ( do_no_rst_count ) {
		if ( !matched ) {
		    (*no_rst_match_count)[list_number]++;
		} else {
		    ++(*rst_count)[list_number];
		}
	    } else {
		++(*rst_count)[list_number];
	    }

	}
    }
    
    if ( no_restr_count != no_restr_count_act) {
	nrerror("Error in Partition_Restraints: bad no_restr_count_act");
    }
#undef NUM_LIST    
}			   



static void Partition_Restraints_ResidueT ( 
    DRestraint *obs_restr, int obs_restr_count, int no_restr_count,
    int **rst_count, Boolean do_no_rst_count, int **no_rst_match_count, 
    Shells *S_T_restr_ass, ID *resID, int num_list )
{

    int		i, no_restr_count_act, atom_idx;
    int		list_number, atom_a_list_number;
    Boolean	double_count, matched;
    ID		id;
    
    *rst_count          = CREATE_NEW( int, num_list );
    *no_rst_match_count = CREATE_NEW( int, num_list );
    
    double_count = OptionTrue( opt_double_count );
    
    no_restr_count_act = 0;
    	
    for (i=0; i<num_list; i++) {
	(*rst_count)[i] = 0;
	(*no_rst_match_count)[i] = 0;
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
        matched = false;

	if ( !obs_restr[i].keep ) {
	    no_restr_count_act++;
	    continue;
	}

	if ( do_no_rst_count ) {
	    if ( S_T_restr_ass[i].S == 0 ) {
	    
		if ( S_T_restr_ass[i].T != 0 ) {	    
		    printf("restraint %d has flags S and T; %d %d\n", i, 
			    S_T_restr_ass[i].S, S_T_restr_ass[i].T);
		    nrerror("Error in Partition_Restraints_ResidueT: bad info ");
		}
		matched = false;
		
	    } else {
	    	matched = true;
	    }
	}
	
	for ( atom_idx = 0; atom_idx < 2; atom_idx++ ) {
	
	    if ( atom_idx == 0 ) {
	        CopyChars( id, obs_restr[i].atom_A.res_nam, size_ID );
	    } else {
	        CopyChars( id, obs_restr[i].atom_B.res_nam, size_ID );
	    }

	    for ( list_number = 0; list_number < num_list; list_number++ ) {
		if ( EQID ( id, resID[ list_number ] ) ) {
		    break;
		}
	    }
	    
	    if ( list_number == num_list ) {
		PRINTIDv( id );
		nrerror ( "Error in Partition_Restraints_ResidueT: wrong residue");
	    }
	    
	    /* Prevent a double count in same list member */

            atom_a_list_number = -1; /* default but just inserted for compiler warnings */
            	    
	    if ( double_count ) {
		if ( atom_idx == 0 ) {
		    atom_a_list_number = list_number;
		} else { 
		    if ( atom_a_list_number == list_number) {
			break;
		    }
		}	    
	    }
	    
	    if ( do_no_rst_count ) {
		if ( !matched ) {
		    (*no_rst_match_count)[list_number]++;
		} else {
		    ++(*rst_count)[list_number];
		}
	    } else {
		++(*rst_count)[list_number];
	    }

	}
    }
    
    if ( no_restr_count != no_restr_count_act) {
	nrerror("Error in Partition_Restraints_ResidueT: bad no_restr_count_act");
    }
}			   


static void Partition_Restraints_Residue ( 
    DRestraint *obs_restr, int obs_restr_count, int no_restr_count,
    int **rst_count, Boolean do_no_rst_count, int **no_rst_match_count, 
    Shells *S_T_restr_ass, Strucset theSet )
{
    
    int		i, no_restr_count_act, atom_idx, num_list;
    int		list_number, atom_a_list_number;
    Boolean	double_count;
    Boolean	matched;
    
    num_list = theSet.molecule->res_count;

    *rst_count          = CREATE_NEW( int, num_list );
    *no_rst_match_count = CREATE_NEW( int, num_list );
    
    double_count = OptionTrue( opt_double_count ); 
    
    no_restr_count_act = 0;
    	
    for (i=0; i<num_list; i++) {
	(*rst_count)[i] = 0;
	(*no_rst_match_count)[i] = 0;
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
        matched = false;
        
	if ( !obs_restr[i].keep ) {
	    no_restr_count_act++;
	    continue;
	}

	if ( do_no_rst_count ) {
	    if ( S_T_restr_ass[i].S == 0 ) {
	    
		if ( S_T_restr_ass[i].T != 0 ) {	    
		    printf("restraint %d has flags S and T; %d %d\n", i, 
			    S_T_restr_ass[i].S, S_T_restr_ass[i].T);
		    nrerror("Error in Partition_Restraints_AtomT: bad info ");
		}
		matched = false;
		
	    } else {
	    	matched = true;
	    }
	}
	
	for ( atom_idx = 0; atom_idx < 2; atom_idx++ ) {
	
	    if ( atom_idx == 0 ) {
		list_number = obs_restr[i].atom_A.res_num_new - 1;
	    } else {
		list_number = obs_restr[i].atom_B.res_num_new - 1;
	    }
	    
	    if ( list_number >= num_list ) {
		printf("%d/n", list_number );
		nrerror ( "Error in Partition_Restraints_ResidueT: wrong residue");
	    }
	    
	    /* Prevent a double count in same list member */
	    
            atom_a_list_number = -1; /* default but just inserted for compiler warnings */

	    if ( double_count ) {
		if ( atom_idx == 0 ) {
		    atom_a_list_number = list_number;
		} else { 
		    if ( atom_a_list_number == list_number) {
			break;
		    }
		}	    
	    }
	    
	    if ( do_no_rst_count ) {
		if ( !matched ) {
		    (*no_rst_match_count)[list_number]++;
		} else {
		    ++(*rst_count)[list_number];
		}
	    } else {
		++(*rst_count)[list_number];
	    }

	}
    }
    
    if ( no_restr_count != no_restr_count_act) {
	nrerror("Error in Partition_Restraints: bad no_restr_count_act");
    }
}


static void Partition_Restraints_Atom ( 
    DRestraint *obs_restr, int obs_restr_count, int no_restr_count,
    int **rst_count, Boolean do_no_rst_count, int **no_rst_match_count, 
    Shells *S_T_restr_ass, int *list_idx,
    int obs_atom_count, 
    Distinfl *dist_infl, int *rev_mapping_table, int dist_count )
{
    
    int		i, j, no_restr_count_act, atom_idx, idx;
    int		list_number, atom_a_list_number;


    Boolean	matched, double_count;
    double_count = OptionTrue( opt_double_count );
    
    *rst_count          = CREATE_NEW( int, obs_atom_count );
    *no_rst_match_count = CREATE_NEW( int, obs_atom_count );
    
    no_restr_count_act = 0;
    	
    for (i=0; i<obs_atom_count; i++) {
	(*rst_count)[i] = 0;
	(*no_rst_match_count)[i] = 0;
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
        matched = false;

	if ( !obs_restr[i].keep ) {
	    no_restr_count_act++;
	    continue;
	}

	if ( do_no_rst_count ) {
	    if ( S_T_restr_ass[i].S == 0 ) {
	    
		if ( S_T_restr_ass[i].T != 0 ) {	    
		    printf("restraint %d has flags S and T; %d %d\n", i, 
			    S_T_restr_ass[i].S, S_T_restr_ass[i].T);
		    nrerror("Error in Partition_Restraints_Atom: bad info ");
		}
		matched = false;
		
	    } else {
	    	matched = true;
	    }
	}
	
	for ( atom_idx = 0; atom_idx < 2; atom_idx++ ) {
	
	    if ( atom_idx == 0 ) {
		idx = obs_restr[i].idx_A;
	    } else {
		idx = obs_restr[i].idx_B;
	    }

	    for ( list_number = 0; list_number < obs_atom_count; list_number++ ) {
	    
		if ( idx == list_idx[list_number] ) {
		    /* printf("Matched atom %d\n", idx); */
		    break;
		}
		
		if ( !do_no_rst_count ) { /* Only look further for the RESTRAINTS */
		    continue;
		}
		
		if ( i == -1 ) {
		    nrerror("Error in Partition_Restraints_Atom -1");
		}
		
		if ( i >= obs_restr_count ) {
		    nrerror("Error in Partition_Restraints_Atom -2");
		}
		
	        j = rev_mapping_table[i];

		if ( j == -1 ) {
		    printf("With %d looking for restraint %d found %d\n", obs_restr_count, i, j);
		    nrerror("Error in Partition_Restraints_Atom 1*");
		}
		
		if ( j >= dist_count ) {
		    nrerror("Error in Partition_Restraints_Atom 2");
		}
		
		/*	
    printf("atom_idx is %d and assignement %d for atom %d\n", atom_idx, dist_infl[j].ass_A, dist_infl[j].A[0].idx );
		*/
		
		if ( atom_idx == 0 ) {
		    if ( dist_infl[j].ass_A == 2) {                /* Pseudo ? */
			if ( dist_infl[j].A[0].idx == list_idx[list_number] ) {
			    break;
			}
		    }
		} else {
		    if ( dist_infl[j].ass_B == 2) {                /* Pseudo ? */
			if ( dist_infl[j].B[0].idx == list_idx[list_number] ) {
			    break;
			}
		    }
		    
		}
	    }
	    
	    if ( list_number == obs_atom_count ) {
		printf("%d %d %d\n", i, atom_idx, idx);
		printf("Doing restraints? %c \n", Boolean_To_Char( do_no_rst_count) );
		nrerror ( "Error in Partition_Restraints_Atom: wrong atom");
	    }
	    
	    /* Prevent a double count in same list member */
	    
            atom_a_list_number = -1; /* default but just inserted for compiler warnings */

	    if ( double_count ) {
		if ( atom_idx == 0 ) {
		    atom_a_list_number = list_number;
		} else { 
		    if ( atom_a_list_number == list_number) {
			break;
		    }
		}	    
	    }
	    
	    if ( do_no_rst_count ) {
		if ( !matched ) {
		    (*no_rst_match_count)[list_number]++;
		} else {
		    ++(*rst_count)[list_number];
		}
	    } else {
		++(*rst_count)[list_number];
	    }

	}
    }
    
    if ( no_restr_count != no_restr_count_act) {
	nrerror("Error in Partition_Restraints: bad no_restr_count_act");
    }
}			   


static void Partition_Restraints_ResAtom ( 
    DRestraint *obs_restr, int obs_restr_count, int no_restr_count,
    int **rst_count, Boolean do_no_rst_count, int **no_rst_match_count, 
    Shells *S_T_restr_ass,
    Strucset theSet, Strucset theSet_ref, 
    Distinfl *dist_infl, int *rev_mapping_table, int dist_count )
{
    
    int		i, j, no_restr_count_act, atom_idx, idx;
    int		list_number, atom_a_list_number, obs_atom_count_ref;
    Boolean	matched, double_count;
    Atom	*atoms, *atoms_ref;
    
    obs_atom_count_ref = theSet_ref.atom_count;
    atoms     = theSet.molecule->atoms;
    atoms_ref = theSet_ref.molecule->atoms;
    
    *rst_count          = CREATE_NEW( int, obs_atom_count_ref );
    *no_rst_match_count = CREATE_NEW( int, obs_atom_count_ref );

    double_count = OptionTrue( opt_double_count );
    
    no_restr_count_act = 0;
    	
    for (i=0; i<obs_atom_count_ref; i++) {
	(*rst_count)[i] = 0;
	(*no_rst_match_count)[i] = 0;
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
        matched = false;

	if ( !obs_restr[i].keep ) {
	    no_restr_count_act++;
	    continue;
	}

	if ( do_no_rst_count ) {
	    if ( S_T_restr_ass[i].S == 0 ) {
	    
		if ( S_T_restr_ass[i].T != 0 ) {	    
		    printf("restraint %d has flags S and T; %d %d\n", i, 
			    S_T_restr_ass[i].S, S_T_restr_ass[i].T);
		    nrerror("Error in Partition_Restraints_ResAtom: bad info ");
		}
		matched = false;
		
	    } else {
	    	matched = true;
	    }
	}
	
	for ( atom_idx = 0; atom_idx < 2; atom_idx++ ) {
	
	    if ( atom_idx == 0 ) {
		idx = obs_restr[i].idx_A;
	    } else {
		idx = obs_restr[i].idx_B;
	    }

	    
	    
	    for ( list_number = 0; list_number < obs_atom_count_ref; list_number++ ) {
	    
		if ( SameResidueT_Atom( atoms[idx], atoms_ref[ list_number ] ) ) {
		    break;
		}
		
		if ( !do_no_rst_count ) { /* Only look further for the RESTRAINTS */
		    continue;
		}
		
		if ( i == -1 ) {
		    nrerror("Error in Partition_Restraints_ResAtom -1");
		}
		
		if ( i >= obs_restr_count ) {
		    nrerror("Error in Partition_Restraints_ResAtom -2");
		}
		
	        j = rev_mapping_table[i];

		if ( j == -1 ) {
		    printf("With %d looking for restraint %d found %d\n", obs_restr_count, i, j);
		    nrerror("Error in Partition_Restraints_ResAtom 1*");
		}
		
		if ( j >= dist_count ) {
		    nrerror("Error in Partition_Restraints_ResAtom 2");
		}
		
		if ( atom_idx == 0 ) {
		    if ( dist_infl[j].ass_A == 2) {                /* Pseudo ? */
			if ( SameResidueT_Atom( atoms[ dist_infl[j].A[0].idx ], 
			     atoms_ref[ list_number ] ) ) {
			    break;
			}
		    }
		} else {
		    if ( dist_infl[j].ass_B == 2) {                /* Pseudo ? */
			if ( SameResidueT_Atom( atoms[ dist_infl[j].B[0].idx ], 
			     atoms_ref[ list_number ] ) ) {
			    break;
			}
		    }
		    
		}
	    }
	    
	    if ( list_number == obs_atom_count_ref ) {
		printf("%d %d %d\n", i, atom_idx, idx);
		printf("Doing restraints? %c \n", Boolean_To_Char( do_no_rst_count) );
		nrerror ( "Error in Partition_Restraints_ResAtom: wrong atom");
	    }
	    
	    /* Prevent a double count in same list member */
	    
            atom_a_list_number = -1; /* default but just inserted for compiler warnings */

	    if ( double_count ) {
		if ( atom_idx == 0 ) {
		    atom_a_list_number = list_number;
		} else { 
		    if ( atom_a_list_number == list_number) {
			break;
		    }
		}	    
	    }
	    
	    if ( do_no_rst_count ) {
		if ( !matched ) {
		    (*no_rst_match_count)[list_number]++;
		} else {
		    (*rst_count)[list_number]++;
		}
	    } else {
		(*rst_count)[list_number]++;
	    }

	}
    }
    
    if ( no_restr_count != no_restr_count_act) {
	nrerror("Error in Partition_Restraints_ResAtom: bad no_restr_count_act");
    }
}			   



static int Categorize_Pairs ( Distval max, int shell_num,
    DRestraint *pair, Distvall *pair_val, int *S_pair_ass, int obs_pair_count,
    int no_pair_count )
    
{
    int i;
    int count, no_pair_count_act;
    Distval avg;
    
    count = 0;
    no_pair_count_act = 0;
    
    avg = GetOptionRVal( opt_avg );
    
    /*
    if ( DEBUG_INFO ) {
	printf("Categorize_Pairs called with %f %d\n", 
	    max, shell_num );
    }
    */
    
    for (i=0; i < obs_pair_count; i++) {
	
	if ( !pair[i].keep  ) {                      /* Skip "deleted" pairs */
	    no_pair_count_act++;
	    continue;
	}
	
	if ( S_pair_ass[i] != 0 ) {        /* Skip pairs already categorized */
	    continue;
	}
	
	if ( Power_Average_Distval_Array ( pair_val[i], avg ) <= max ) {
	    S_pair_ass[i] = shell_num;
	    count++;
	} 
    }
        
    if ( no_pair_count_act != no_pair_count ) {
	printf("%d %d \n", no_pair_count_act, no_pair_count);
	nrerror("Error in Categorize_Pairs: bad no_pair_count_act");
    }
    
    return ( count );
}	    




static void Fill_Int_Array( int *m, int value, int nrl, int nrh)
{
    int r;
    
    for (r=nrl; r<=nrh; r++) {	
	m[r] = value;
    }	
    
    return;
}

static void Fill_Boolean_Array( Boolean *m,  Boolean value, int nrl, int nrh)
{
    int r;
    
    for (r=nrl; r<=nrh; r++) {	
	m[r] = value;
    }	
    
    return;
}

static void Fill_2Int_Array( Shells *m, int value1, int value2, 
				    int nrl, int nrh)
{
    int r;
    
    for (r=nrl; r<=nrh; r++) {	
	m[r].S = value1;
	m[r].T = value2;
    }	
    
    return;
}




static int Match_Rsts_Pairs ( int s,
    DRestraint *obs_restr, Shells *S_T_restr_ass, int obs_restr_count, int no_restr_count,
    DRestraint *pair, int *S_pair_ass, int obs_pair_count, 
    Distinfl *dist_infl, int *rev_mapping_table, 
    Boolean *pair_matched ) 
{
#define SMALL_FLOAT 0.0001

    int r, t, i, j;
    int smaller_count, larger_count;
    int shell_count, no_restr_count_act, bin_offset;
    int number_shells_rst, number_shells_prs;
    Distval minimum_distance_rst, maximum_distance_rst;
    Boolean found;
    Distval bin_size;

    number_shells_rst    = GetOptionIVal( opt_numb_shells_rst );
    minimum_distance_rst = GetOptionRVal( opt_min_dist_rst );
    maximum_distance_rst = GetOptionRVal( opt_max_dist_rst );
    
    number_shells_prs    = GetOptionIVal( opt_numb_shells_prs );
        
    shell_count        = 0;
    no_restr_count_act = 0;
    smaller_count      = 0;
    larger_count       = 0;

    bin_size = (maximum_distance_rst - minimum_distance_rst) / number_shells_rst;
    bin_offset =  (int) ( (minimum_distance_rst / bin_size) + SMALL_FLOAT - 1 );

    if ( terseness < 6 ) {
        printf("Matching shell %d\n", s);
    }
    
    if ( ( s < 1 ) || ( s > number_shells_prs ) ) {
	printf("given shell is %d\n", s);
	nrerror("Error in Match_Rst_Pair: bad shell info ");
    }
    
    for (i=0; i<obs_restr_count; i++) {
    
	if ( !obs_restr[i].keep  ) {         /* Skip "deleted" restraints */
	    no_restr_count_act++;
	    continue;
	}
	
	if ( S_T_restr_ass[i].S != 0 ) {     /* Skip already categorized */
	    continue;
	}
	
	if ( S_T_restr_ass[i].T != 0 ) {	    
	    printf("restraint %d has flags S and T; %d %d\n", i, 
	    S_T_restr_ass[i].S, S_T_restr_ass[i].T);
	    nrerror("Error in Match_Rst_Pair: bad info ");
	}
	
	found = false;
	
	r = rev_mapping_table[i];
	
	for (j=0; j< obs_pair_count; j++) {
	
	    if ( !pair[j].keep  ) {                  /* Skip "deleted" pairs */
		continue;
	    }

	    if ( S_pair_ass[j] != s  ) {     /* Skip pairs in other category */
		continue;
	    }
	    
	    if ( pair_matched[j] ) {          /* Skip pairs already matched */
		continue;
	    }
	    
	    if ( Same_Partners_DRestraint( pair[j], obs_restr[i], 
			    dist_infl, r ) ) {
		found = true;
		pair_matched[j] = true;
		
		/*
		printf("Match P%5d - R%5d ", j, i);
		printf("%3d ", pair[j].atom_A.res_num_orig);
		PRINTID( pair[j].atom_A.res_nam  ); printf(" ");		 	
		PRINTID( pair[j].atom_A.atom_int ); printf(" ");	 	
		printf("%3d ", pair[j].atom_B.res_num_orig);
		PRINTID( pair[j].atom_B.res_nam  ); printf(" ");		 	
		PRINTID( pair[j].atom_B.atom_int ); printf("        ");	 	
    
		printf("%3d ", obs_restr[i].atom_A.res_num_orig);
		PRINTID( obs_restr[i].atom_A.res_nam  ); printf(" ");		 	
		PRINTID( obs_restr[i].atom_A.atom_int ); printf(" ");	 	
		printf("%3d ", obs_restr[i].atom_B.res_num_orig);
		PRINTID( obs_restr[i].atom_B.res_nam  ); printf(" ");		 	
		PRINTID( obs_restr[i].atom_B.atom_int ); printf("\n");
		*/
		
		break;
	    }
	}
	
	if ( found ) {
	    
	    S_T_restr_ass[i].S = s;
	    
	    t = (int) ( (obs_restr[i].up_bnd / bin_size) - bin_offset + SMALL_FLOAT );
	    
	    /*
	    printf("up_bnd %8.3f bin_size %8.3f -> t %d\n", obs_restr[i].up_bnd, 
		    bin_size, t);
	    */
	        
	    if ( t < 1 ) {
		t = 1;
		smaller_count++;
	    }
	    
	    if ( t > number_shells_rst ) {
		t = number_shells_rst;
		larger_count++;
	    }
	    	    
	    S_T_restr_ass[i].T = t;
	    	    	
	    shell_count++;
	}
    }


    if ( no_restr_count_act != no_restr_count ) {
	printf("%d %d \n", no_restr_count_act, no_restr_count );
	nrerror("Error in Categorize_Pairs: bad no_restr_count_act");
    }


    if ( ( smaller_count > 0 ) && ( terseness < 6 ) ) {
	printf("Found %d restraints smaller than minimum distance\n",smaller_count);
    }
    
    if ( (  larger_count > 0 ) && ( terseness < 6 ) ){
	printf("Found %d restraints larger than maximum distance\n",  larger_count);
    }
    
    /*
    printf("Total of %d restraints classified\n", shell_count);
    */
#undef SMALL_FLOAT
    
    return ( shell_count );
}



static Boolean Same_Partners_DRestraint ( DRestraint P, DRestraint Q, 
Distinfl *dist_infl, int r )
{

/* Checks if atom A of Drestraint P and Q are equal and
 *           atom B of Drestraint P and Q are equal (or vica versa)
 * P is theoretical and Q experimental so
 * If atom A in Q is expandable to atom A in P it's a match (or vica versa)
 */
    Boolean match_1, match_2;
    
    /*  Early return condition, not worthwhile nop
    if ( ( ( P.idx_A == Q.idx_A ) && ( P.idx_B == Q.idx_B ) ) ||
	 ( ( P.idx_A == Q.idx_B ) && ( P.idx_B == Q.idx_A ) )  ) {
	return ( true );
    }

    if ( r == -1 ) {
	nrerror("Error in Same_Partners_DRestraint 1");
    }
    if ( r >= dist_count ) {
	nrerror("Error in Same_Partners_DRestraint 2");
    }
    */
    
    if ( P.idx_A == Q.idx_A ) {
	match_1 = true;
    } else {
	if ( dist_infl[r].ass_A == 2) {                /* Pseudo ? */
	    if ( ( dist_infl[r].A[0].idx == P.idx_A ) ||
		 ( dist_infl[r].A[1].idx == P.idx_A ) ) {
		match_1 = true;
	    } else {
		match_1 = false;
	    }
	} else {
	    match_1 = false;
	}
    }

    if ( P.idx_B == Q.idx_B ) {
	match_2 = true;
    } else {
	if ( dist_infl[r].ass_B == 2) {                /* Pseudo ? */
	    if ( ( dist_infl[r].B[0].idx == P.idx_B ) ||
		 ( dist_infl[r].B[1].idx == P.idx_B ) ) {
		match_2 = true;
	    } else {
		match_2 = false;
	    }
	} else {
	    match_2 = false;
	}
    }


    if ( match_1 && match_2 ) {
	return ( true );
    }



    if ( P.idx_B == Q.idx_A ) {
	match_1 = true;
    } else {
	if ( dist_infl[r].ass_A == 2 ) {                /* Pseudo ? */
	    if ( ( dist_infl[r].A[0].idx == P.idx_B ) ||
		 ( dist_infl[r].A[1].idx == P.idx_B ) ) {
		match_1 = true;
	    } else {
		match_1 = false;
	    }
	} else {
	    match_1 = false;
	}
    }


    if ( P.idx_A == Q.idx_B ) {
	match_2 = true;
    } else {
	if ( dist_infl[r].ass_B == 2) {                /* Pseudo ? */
	    if ( ( dist_infl[r].B[0].idx == P.idx_A ) ||
		 ( dist_infl[r].B[1].idx == P.idx_A ) ) {
		match_2 = true;
	    } else {
		match_2 = false;
	    }
	} else {
	    match_2 = false;
	}
    }

    if ( match_1 && match_2 ) {
	return ( true );
    } else {
	return ( false );	
    }
}
    


static void ShowSummary (   
		    int obs_restr_count,   int obs_pair_count, 
		    Shells *S_T_restr_ass, int *S_pair_ass)
{
    int s, t;
    int observed, expected, completeness, completeness_total;
    int observed_total, expected_total;
    Distval bin_rst, low_rst, high_rst;
    Distval bin_prs, low_prs, high_prs;
    int number_shells_rst, number_shells_prs;
    Distval minimum_distance_rst, maximum_distance_rst;
    Distval minimum_distance_prs, maximum_distance_prs;
    
    observed_total = 0;
    expected_total = 0;

    number_shells_rst    = GetOptionIVal( opt_numb_shells_rst );
    minimum_distance_rst = GetOptionRVal( opt_min_dist_rst );
    maximum_distance_rst = GetOptionRVal( opt_max_dist_rst );
    
    number_shells_prs    = GetOptionIVal( opt_numb_shells_prs );
    minimum_distance_prs = GetOptionRVal( opt_min_dist_prs );
    maximum_distance_prs = GetOptionRVal( opt_max_dist_prs );
    
    bin_rst = (maximum_distance_rst - minimum_distance_rst) / number_shells_rst;
    bin_prs = (maximum_distance_prs - minimum_distance_prs) / number_shells_prs;
    
    printf("\nSUMMARY\n");
    
    printf("Shell   Expect ");
    for (t=1; t<= number_shells_rst; t++) {

	low_rst  = minimum_distance_rst + (t-1)*bin_rst;
	printf("%3.1f- ", low_rst);
    }
    printf(" TOTAL SAT CUM\n");

    printf("                ");
    for (t=1; t<= number_shells_rst; t++) {

	high_rst  = minimum_distance_rst + t*bin_rst;
	printf("%3.1f  ", high_rst);
    }
    printf("\n");
   
    for (s=1; s<= number_shells_prs; s++) {

	low_prs  = minimum_distance_prs + (s-1)*bin_prs;
	high_prs = low_prs + bin_prs;
	printf("%3.1f-%3.1f ", low_prs, high_prs);
	
	expected = Count_Pair_S (s, S_pair_ass, obs_pair_count);
	printf("%6d ", expected );
	expected_total += expected;
	
	for (t=1; t<= number_shells_rst; t++) {
	    printf("%4d ", Count_Restraint_S_T ( s, t, 
					    S_T_restr_ass, obs_restr_count) );
	}
	
	observed = Count_Restraint_S (s, S_T_restr_ass, obs_restr_count);
	observed_total += observed;
	

	if ( expected != 0 ) {
	    completeness = (int) (100.0*observed / expected + 0.5);    
	} else {
	    completeness = (int) NOT_A_NUMBER;
	}


	if ( expected_total != 0 ) {
	    completeness_total = (int) (100.0*observed_total / expected_total + 0.5);    
	} else {
	    completeness_total = (int) NOT_A_NUMBER;
	}

	printf("%6d %3d %3d\n", observed, completeness, completeness_total );
    }
}		    



static void Show_Summary_AQUA ( FILE *fp,   
    int obs_restr_count,   int obs_pair_count, 
    Shells *S_T_restr_ass, int *S_pair_ass)
    
{
    int s, t;
    int observed, expected, completeness, completeness_total;
    int observed_total, expected_total;
    Distval bin_rst, low_rst;
    Distval bin_prs, low_prs, high_prs;
    int number_shells_rst, number_shells_prs;
    Distval minimum_distance_rst, maximum_distance_rst;
    Distval minimum_distance_prs, maximum_distance_prs;
    
    observed_total = 0;
    expected_total = 0;

    number_shells_rst    = GetOptionIVal( opt_numb_shells_rst );
    minimum_distance_rst = GetOptionRVal( opt_min_dist_rst );
    maximum_distance_rst = GetOptionRVal( opt_max_dist_rst );
    
    number_shells_prs    = GetOptionIVal( opt_numb_shells_prs );
    minimum_distance_prs = GetOptionRVal( opt_min_dist_prs );
    maximum_distance_prs = GetOptionRVal( opt_max_dist_prs );
    
    bin_rst = (maximum_distance_rst - minimum_distance_rst) / number_shells_rst;
    bin_prs = (maximum_distance_prs - minimum_distance_prs) / number_shells_prs;
    
    fprintf( fp, "\n" );
    fprintf( fp, "$ DATA  COMPLETENESS_PER_SHELL\n" );
    fprintf( fp, "$ TABLE #RECORDS %d #ITEMS %d\n", 
		 number_shells_prs, 6 + number_shells_rst);
    
    fprintf( fp, "  LowP   UppP     NP " );
    for (t=1; t<= number_shells_rst; t++) {
	low_rst  = minimum_distance_rst + (t-1)*bin_rst;
	fprintf( fp, "%3.1f- ", low_rst );
    }
    fprintf( fp, "  N_RM      C  C_cum\n" );

    for (s=1; s<= number_shells_prs; s++) {

	low_prs  = minimum_distance_prs + (s-1)*bin_prs;
	high_prs = low_prs + bin_prs;
	fprintf( fp, "%6.2f %6.2f ", low_prs, high_prs );
	
	expected = Count_Pair_S (s, S_pair_ass, obs_pair_count);
	fprintf( fp, "%6d ", expected );
	expected_total += expected;
	
	for (t=1; t<= number_shells_rst; t++) {
	    fprintf( fp, "%4d ", Count_Restraint_S_T ( s, t, 
					    S_T_restr_ass, obs_restr_count) );
	}
	
	observed = Count_Restraint_S (s, S_T_restr_ass, obs_restr_count);
	observed_total += observed;
	
	if ( expected != 0 ) {
	    completeness = (int) (100.0*observed / expected + 0.5);    
	} else {
	    completeness = (int) NOT_A_NUMBER;
	}

	if ( expected_total != 0 ) {
	    completeness_total = (int) (100.0*observed_total / expected_total + 0.5);    
	} else {
	    completeness_total = (int) NOT_A_NUMBER;
	}

	fprintf( fp, "%6d %6d %6d\n", observed, completeness, completeness_total );
    }
    fprintf( fp, "$ END\n" );
}		    



static int Count_Restraint_S ( int shells_s, Shells *S_T_restr_ass, 
				int obs_restr_count )
{
    int i, count;
    
    count = 0;
    
    for (i=0; i<= obs_restr_count; i++) {
	if ( S_T_restr_ass[i].S == shells_s ) {
	    count++;
	}
    }
    
    return ( count );
}


static int Count_Restraint_S_T ( int shells_s, int shells_t, 
			Shells *S_T_restr_ass, int obs_restr_count )
{
    int i, count;
    
    count = 0;
    
    for (i=0; i<= obs_restr_count; i++) {
	if ( ( S_T_restr_ass[i].S == shells_s ) &&
	     ( S_T_restr_ass[i].T == shells_t ) ) {
	    count++;
	}
    }
    
    return ( count );
}


static int Count_Pair_S ( int shells_s, int *S_pair_ass, int obs_pair_count )
{
    int i, count;
    
    count = 0;
    
    for (i=0; i<= obs_pair_count; i++) {
	if ( S_pair_ass[i] == shells_s ) {
	    count++;
	}
    }
    
    return ( count );
}


static char Boolean_To_Char(Boolean boo)
{
    
    if (boo) 
	return('T'); 
    else 
	return('F');
}



static Boolean Same_Residue (Atom a, Atom b)
{
    if ( NOTEQID(a.chain, b.chain) ) {
	return ( false );
    }
    
    if ( a.res_num_orig == b.res_num_orig ) {
	return ( true );
    } else {
	return ( false );
    }
}



static void Match_SSA ( Strucset *setptr, Distinfl *dist_infl, int dist_count  )  
{
/* match triplets present in the mol with ssa restraints and
    write results to temporary file  */

    int	    i, j, k, count, idx, ass;
    int	    stereo_count, total_ssa;
    Psatom  *stereo_list;
    Boolean *stereo_assigned;

    FILE    *fp;
    char    File[]="AquaCompl.temp";
    Mol *molptr;
    
    molptr = setptr->molecule;
    total_ssa = 0;
    
    /* generate the required "stereo" information about equivalent atoms */

    if ( OptionTrue( opt_sterlib ) )
	ModOption( opt_rstinp, "user" );  
	/* force user libraries; cf. AquaPseudo */

    stereo_list = ProcessStereoLIB( &stereo_count, molptr, "StereoLIB-aqua"  );
    if ( !stereo_list || stereo_count == 0 )
    {
	puts( " * No stereo definitions!" );
    }

    stereo_assigned = CREATE_NEW( Boolean, stereo_count );
    
    for (i=0; i< stereo_count; i++) {
	
	stereo_assigned[i] = false;    
    
	for (j=0; j< dist_count; j++) {
	
	    for (k=0; k<2; k++) {
		
		if ( k == 0 ) {
		    count = dist_infl[j].count_A;
		    ass   = dist_infl[j].ass_A;
		} else {
		    count = dist_infl[j].count_B;
		    ass   = dist_infl[j].ass_B;
		}

		if ( ( count == 1 ) || ( ass == 2 ) ) {
		    continue;
		}

		if ( k == 0 ) {
		    idx   = dist_infl[j].A[ 0 ].idx;
		} else {
		    idx   = dist_infl[j].B[ 0 ].idx;
		}

		if ( Same_Atom ( molptr->atoms[ idx ], stereo_list[i] ) ) {
		
		    if ( stereo_assigned[i] ) {
			nrerror("Error in Match_SSA: already assigned");
		    }
		    
		    total_ssa++;
		    stereo_assigned[i] = true;
		    break;
		}
		
	    }
	    
	    if ( stereo_assigned[i] ) {
		break;
	    }
	}
    }
    
    if ( ( fp = fopen( File, "w" ) ) == NULL ) {
        printf( " * File <%s> could not be opened in mode <%s>\n", File, "w" );
	nrerror("Error in Match_SSA");
    }

    if ( terseness < 6 ) {
        printf("Opened temporary file %s\n", File);
    }
	    
    Show_Type_SSA_info_AQUA ( fp,
	    stereo_list, stereo_assigned, stereo_count, total_ssa );
    
    if ( OptionTrue( opt_detail ) ) {
	Show_SSA_info_AQUA ( fp, stereo_list, stereo_assigned, stereo_count );
    }
    
    fclose ( fp );
}



static Boolean Same_Atom ( Atom atom_A, Psatom psatom_B ) 
{    

/* Checks if atom_A is the first atom of the form atoms psatom_B */

    if ( atom_A.res_num_orig != psatom_B.res_num_orig ) {
	return ( false );
    }

    if ( NOTEQID(atom_A.atom_int, psatom_B.form_nam[ 0 ]) ) {
	return ( false );
    } else {
	if ( EQID(atom_A.chain, psatom_B.chain ) ) {
	    return ( true );
	} else {
	    return ( false );
	}
    }  
}


static void Show_Type_SSA_info_AQUA ( FILE *fp,
	Psatom *stereo_list, 
	Boolean *stereo_assigned, int stereo_count, int total_ssa )

{
#define NUM_TYPES 6
    int	    i, j;
    int	    *theor, *exper;
    int	    total_theor, total_exper;
    Distval   completeness;
    
    total_theor = 0;
    total_exper = 0;
    
    theor = CREATE_NEW( int, NUM_TYPES+1 );
    exper = CREATE_NEW( int, NUM_TYPES+1 );
    
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  SSA_PER_TYPE\n" );
    fprintf( fp, "$ REMARK   pseudoatom type is set according to the following table:\n" );
    fprintf( fp, "$ REMARK    The number between brackets is the number of atoms within.\n" );
    fprintf( fp, "$ REMARK  		0 undefined  (*)\n" );
    fprintf( fp, "$ REMARK  		1 CH2 or NH2 (2) 		2 methyl or NH3 (3)\n" );
    fprintf( fp, "$ REMARK  		3 double CH2/NH2 (4)		4 double methyl (6)\n" );
    fprintf( fp, "$ REMARK  		5 Aromat with 2 H (2) 		6 Aromat with 4 H (4)\n" );
    fprintf( fp, "$ REMARK  In this file pseudo's of type 0 and 2 can not be used\n" );
    fprintf( fp, "$ REMARK  Type 6 depends on type 5 and is not (yet) properly scored\n" );
    fprintf( fp, "$ REMARK  Theor : stereospecific groups present\n" );
    fprintf( fp, "$ REMARK  Exper : stereospecific assigned according to restraints\n" );
    fprintf( fp, "$ REMARK  Perc  : percentage of the Exper to the Theor\n" );
    fprintf( fp, "$ TABLE #RECORDS 7 #ITEMS 4\n" );
    fprintf( fp, "Type    Theor  Exper   Perc\n" );
    
    for (i=1; i<= NUM_TYPES; i++) {
    
	theor[i] = 0;

	for (j=0; j< stereo_count; j++) {
	    if ( stereo_list[j].pseudo_typ == i ) {
		theor[i]++;
		if ( stereo_assigned[j] ) {
		    exper[i]++;
		}
	    }
	}
	
	if ( theor[i] == 0 ) {
	    completeness = NOT_A_NUMBER;
	} else {
	    completeness = 100.0 * exper[i] / theor[i];
	}	
	
	fprintf( fp, "%6d %6d %6d %6.1f\n", i, theor[i], exper[i], 
		completeness );
	    
	total_theor += theor[i];
	total_exper += exper[i];
    }

    if ( total_ssa != total_exper ) {
	nrerror("Error in Show_Type_SSA_info_AQUA");
    }
    
    if ( total_theor == 0 ) {
	completeness = NOT_A_NUMBER;
    } else {
	completeness = 100.0 * total_exper / total_theor;
    }	
    
    fprintf( fp, "%6s %6d %6d %6.1f\n", "All",
		total_theor, total_exper, completeness );
    
    fprintf( fp, "$ END\n" );

#undef NUM_TYPES
}


static void Show_SSA_info_AQUA ( FILE *fp,
	    	Psatom *stereo_list, 
		Boolean *stereo_assigned, int stereo_count )
{
    int	    i;
        
    fputc( '\n', fp );
    fprintf( fp, "$ DATA  SSA_PER_ATOM\n" );
    fprintf( fp, "$ REMARK  PS    : Pseudo atom representing the stereospecific group\n" );
    fprintf( fp, "$ REMARK  AtomA : First Atom of the stereospecific group\n" );
    fprintf( fp, "$ REMARK  SSA   : Assignement of the group; assigned if true (T)\n" );
    fprintf( fp, "$ TABLE #RECORDS %d #ITEMS 8\n", stereo_count );
    fprintf( fp, "Chn Resnm Res  PS Type AtomA AtomB SSA\n" );
    
    for (i=0; i< stereo_count; i++) {
    
        WriteFPsatom ( fp, stereo_list[i] );      /* 7 items */
	
	fprintf( fp, "%c\n", Boolean_To_Char( stereo_assigned[i] ) ); 	    
    }    
    fprintf( fp, "$ END\n" );
}




static Boolean Concatenate_file_A_To_B ( char *file_a,  char *file_b )
{
/* Append file A to file B and remove the file A */

    char *cmmnd;
    
    cmmnd = CREATE_NEW( char, 256 );


    sprintf( cmmnd, "cat %s %s > temp", file_b, file_a );
    if ( terseness < 6 ) {
        printf ("Executing csh command: %s\n", cmmnd);
    }
    if ( system( cmmnd ) != 0 ) {
	return ( false );
    }

    sprintf( cmmnd, "\\mv -f temp %s", file_b );
    if ( terseness < 6 ) {
        printf ("Executing csh command: %s\n", cmmnd);
    }
    if ( system( cmmnd ) != 0 ) {
	return ( false );
    }
    
    sprintf( cmmnd, "\\rm -f %s", file_a );
    if ( terseness < 6 ) {
        printf ("Executing csh command: %s\n", cmmnd);
    }
    if ( system( cmmnd ) != 0 ) {
	return ( false );
    }

    return ( true );
}



static void WriteFPsatom ( FILE *fp, Psatom psatom ) 
{
    size_t  l;

    fputc( CHAIN, fp );
    if ( ( l = LenChars( psatom.chain, size_ID ) ) )
    {
	FPrintChars( fp, psatom.chain, l );
        fputc( ' ', fp );
    }
    else
        fprintf( fp, "  " );

    FPRINTID( fp, psatom.res_nam );
    fprintf( fp, " %3i ", psatom.res_num_orig );

    fprintf( fp, "  " );

    FPRINTID( fp, psatom.pseudo_nam );
    
    fprintf( fp, " %1i ", psatom.pseudo_typ );
    
    FPRINTID( fp, psatom.form_nam[0] );
    fprintf( fp, " " );
    FPRINTID( fp, psatom.form_nam[1] );
  
}


static Boolean SameResidueT_Atom ( Atom A, Atom B )
{
/* Checks if the two atoms have the same atom and residue name */

    if ( EQID ( A.atom_int, B.atom_int ) &&
	 EQID ( A.res_nam,  B.res_nam ) ) {
	return ( true );
    } else {
	return ( false );	
    }
}

static Boolean Check_Options ( void )
{

#define SMALL_FLOAT 0.0001
#define MAX_SHELLS   99
#define MAX_DISTANCE 999.9

    Boolean check;
    
    check = true;
    
    if ( ( GetOptionIVal( opt_numb_shells_prs ) < 1 ) || 
         ( GetOptionIVal( opt_numb_shells_prs ) >= MAX_SHELLS ) ) {
	check = false;
	printf("opt_numb_shells_prs %d\n",GetOptionIVal(opt_numb_shells_prs));
    }

    if ( ( GetOptionIVal( opt_numb_shells_rst ) < 1 ) || 
         ( GetOptionIVal( opt_numb_shells_rst ) >= MAX_SHELLS ) ) {
	check = false;
	printf("opt_numb_shells_rst %d\n",GetOptionIVal(opt_numb_shells_rst));
    }

    if ( GetOptionRVal( opt_min_dist_prs ) <= SMALL_FLOAT ) {
	check = false;
	printf("opt_min_dist_prs %8.3f\n",GetOptionRVal( opt_min_dist_prs ));
    }

    if ( GetOptionRVal( opt_max_dist_prs ) >= MAX_DISTANCE ) {
	check = false;
	printf("opt_min_dist_prs %8.3f\n",GetOptionRVal( opt_min_dist_prs ));
	printf("opt_max_dist_prs %8.3f\n",GetOptionRVal( opt_max_dist_prs ));
    }

    if ( GetOptionRVal( opt_max_dist_prs ) <= GetOptionRVal( opt_min_dist_prs ) ) {
	check = false;
	printf("opt_min_dist_prs %8.3f\n",GetOptionRVal( opt_min_dist_prs ));
	printf("opt_max_dist_prs %8.3f\n",GetOptionRVal( opt_max_dist_prs ));
    }

    if ( GetOptionRVal( opt_min_dist_rst ) <= SMALL_FLOAT ) {
	check = false;
	printf("opt_min_dist_rst %8.3f\n",GetOptionRVal( opt_min_dist_rst ));
    }

    if ( GetOptionRVal( opt_max_dist_rst ) >= MAX_DISTANCE ) {
	check = false;
	printf("opt_min_dist_rst %8.3f\n",GetOptionRVal( opt_min_dist_rst ));
	printf("opt_max_dist_rst %8.3f\n",GetOptionRVal( opt_max_dist_rst ));
    }

    if ( GetOptionRVal( opt_max_dist_rst ) <= GetOptionRVal( opt_min_dist_rst ) ) {
	check = false;
	printf("opt_min_dist_rst %8.3f\n",GetOptionRVal( opt_min_dist_rst ));
	printf("opt_max_dist_rst %8.3f\n",GetOptionRVal( opt_max_dist_rst ));
    }

    if ( ( GetOptionRVal( opt_avg ) >  2.0 ) || 
         ( GetOptionRVal( opt_avg ) < -9.0 ) ) {
	check = false;
	printf("opt_avg %8.3f\n",GetOptionRVal( opt_avg ));
    }

    if ( strlen( GetOptionCVal(opt_atominp) ) == 0 ) {
	check = false;
	printf("opt_atominp is empty\n");	     
    }

    if ( GetOptionRVal( opt_thresh ) < 0.0 ) {
	printf("WARNING in Check_Options: opt_thresh %8.3f\n",GetOptionRVal( opt_thresh ));
    }
    if ( terseness < 9 ) {
        printf ("Checked all options\n");
    }
	    
    return ( check );
    
#undef SMALL_FLOAT
#undef MAX_SHELLS
#undef MAX_DISTANCE
}
